NotebookDirectory[]

04-05 Kriging

({{Copyright © A. I. McLeod, 2015, }})

Conditional Mean and Variance in Normal Distribution

Let u and v be jointly normally distributed multivariate random vectors of dimensions (Subscript[m, u]) and (Subscript[m, v]) with means (Subscript[μ, u]) and (Subscript[μ, v]), covariance matricies (Subscript[C, u,u]) and (Subscript[C, v,v]) and cross-covariances, (Subscript[C, u,v])=𝔈{(u-(Subscript[μ, u]))((v-(Subscript[μ, v])))'}=(Subscript[C, v,u]^′).

(({{u}, {v}})) ~ Gau((({{(Subscript[μ, u])}, {(Subscript[μ, v])}})),(({{(Subscript[C, u,u]), (Subscript[C, u,v])}, {(Subscript[C, v,u]), (Subscript[C, v,v])}})))

Textbooks that discuss the multivariate normal distribution derive the conditional mean and variance in this situation.  See also, http://en.wikipedia.org/wiki/Multivariate_normal_distribution   LINK . So we have,

𝔈{u|v=(Subscript[V, o])}=(Subscript[μ, u])+(Subscript[C, u,v])(Subscript[C, v,v]^(-1))((Subscript[V, o])-(Subscript[μ, v]))

and

Var(u|v=(Subscript[V, 0]))=(Subscript[C, u,u])-(Subscript[C, u,v])(Subscript[C, v,v]^(-1))(Subscript[C, v,u])

Eqn. () is the multivariate regression of u on v with prediction variance given by ().

(Subscript[C, u,v])(Subscript[C, v,v]^(-1)) are the regression coefficients for the regression of U on V.

These results are easily derived using the Sherman-Morrison-Woodbury result for the inverse of a partitioned matrix.

Spline Regression

Spline regression is a good way to estimate and remove spatial trends from the data when developing BHM models.

Some Univariate Basis Functions

The constant term is not shown.  Note in R, the function poly() constructs orthogonal polynomials - these are computed using the Gram-Schmidt orthogonalization method.

#Source: BasesFunctionsSplines
#Bases functions: splines & polynomials
library(splines)
library(lattice)
x<-seq(from=1, to=2, length=100)
n<-length(x)
B<-c(bs(x, 3))
N<-c(ns(x, 3))
P<-c(poly(x,3,raw=TRUE))
O<-c(poly(x,3))
y<-c(P,O,B,N)
Bnames<-rep("B-splines", 3*n)
Nnames<-rep("Natural-splines", 3*n)
Pnames<-rep("Polynomial", 3*n)
Onames<-rep("Orthogonal..", 3*n)
type<-ordered(c(Pnames,Onames,Bnames,Nnames),levels=c("Polynomial","Orthogonal..","B-splines","Natural-splines"))
degree<-ordered(rep(rep(1:3, rep(n,3)),4))
DF<-data.frame(x=x, y=y, type=type, degree=degree)
xyplot(y~x|type*degree, data=DF, type="l",
    scales=list(y=list(relation="free", draw=FALSE)),
    main="Basis functions",lwd=3
    )

GAG Urine data example

In the script below we fit best AIC and BIC spline regression models. The AIC model fit is not shown here but it has df=22 and is only slightly different from BIC.

#source: SplineRegressionGAGAICBIC.R
#Regression using natural splines with AIC/BIC
#
require("splines")
require("MASS")
data(GAGurine)
names(GAGurine)
#Natural Spline Regression, select with AIC/BIC
#note build-in lm methods AIC(), BIC(), logLik()
dfs <- 10:30
numdfs <- length(dfs)
bics <- aics <- numeric(numdfs)
names(aics) <- dfs
for (i in 1:numdfs) {
  ans <- lm(GAG~ns(Age, df=dfs[i]),data=GAGurine)
  aics[i] <- AIC(ans)
  bics[i] <- BIC(ans)
}
#best AIC model
dfAIC <- dfs[which.min(aics)]
ans<-lm(GAG~ns(Age, df=dfAIC),data=GAGurine)
aicans <- round(AIC(ans),2)
x<-seq(from=min(GAGurine$Age), to=max(GAGurine$Age), length=100)
y<-predict(ans, newdata=data.frame(Age=x))
plot(GAGurine$Age, GAGurine$GAG, xlab="Age", ylab="GAG")
lines(x,y, lwd=2, col="red")
TitleA <- paste("Natural Spline Regression, df =", dfAIC)
TitleB <- paste(TitleA, ", AIC =")
titleMain <- paste(TitleB, aicans)
title(main=titleMain)
#
#best BIC model
dfBIC <- dfs[which.min(bics)]
ans<-lm(GAG~ns(Age, df=dfBIC),data=GAGurine)
bicans <- round(BIC(ans),2)
x<-seq(from=min(GAGurine$Age), to=max(GAGurine$Age), length=100)
y<-predict(ans, newdata=data.frame(Age=x))
plot(GAGurine$Age, GAGurine$GAG, xlab="Age", ylab="GAG")
lines(x,y, lwd=2, col="red")
TitleA <- paste("Natural Spline Regression, df =", dfBIC)
TitleB <- paste(TitleA, ", BIC =")
titleMain <- paste(TitleB, bicans)
title(main=titleMain)

Tensor Splines

Consider the case of two inputs (Subscript[X, 1]) and (Subscript[X, 2]). Denote a set of basis functions for (Subscript[X, 1]) and (Subscript[X, 2]) by (Subscript[h, 1,j])((Subscript[X, 1])),j=1,…,(Subscript[M, 1]) and (Subscript[h, 2,k])((Subscript[X, 2])), k=1,…,(Subscript[M, 2]).  Then the tensor product basis for X=((Subscript[X, 1]),(Subscript[X, 2])) is defined by (Subscript[g, j,k])(X)=(Subscript[h, 1,j])((Subscript[X, 1]))(Subscript[h, 2,k])((Subscript[X, 2])).  For B-splines the tensor basis is illustrated below.  Notice the univariate basis functions for (Subscript[X, 1]) and (Subscript[X, 2]) may have different shapes since they may be defined on different intervals.

"Figure 5.10: A tensor product basis of B-splines, showing some selected pairs. Each two-dimensional function is the tensor product of the corresponding one dimensional marginals."  HTF

Comparison: Additive splines vs. tensor splines for mixture dataset

The additive splines make the strong assumption of additivity and so lose some flexibility if there is real mutual dependence or interaction effects between the variables.  With additive splines the dependence component in the model is described by,

g(X)=(∑_(j=1)^(Subscript[M, 1]))(Subscript[θ, j])(Subscript[g, j])((Subscript[X, 1]))+(∑_(k=1)^(Subscript[M, 2]))(Subscript[θ, k])(Subscript[g, k])((Subscript[X, 2]))

In this case there would be (Subscript[M, 1])+(Subscript[M, 2]) parameters.  If constant terms are present, then the degrees of freedom are (Subscript[M, 1])+(Subscript[M, 2])-1.

"Figure 5.11: The simulation example of Figure 2.1. This panel shows the decision boundary of an additive logistic regression model, using natural splines in each of the two coordinates (total df = 1+(4 -1) +(41) = 7).  The broken purple boundary is the Bayes decision boundary for this problem."   HTF

Tensor product splines provide a much more flexible fit but this doesn't usually mean better predictions.  With tensor product splines our model component can be written,

g(X)=(∑_(j=1)^(Subscript[M, 1]))(∑_(k=1)^(Subscript[M, 2]))(Subscript[θ, j,k])(Subscript[g, j,k])(X)

and there are (Subscript[M, 1])×(Subscript[M, 2]) parameters.  For more than 2 inputs, this approach rapidly becomes infeasible due to the rapid increase in parameters.

"Figure 5.11: The simulation example of Figure 2.1. This panel shows the results of using a tensor product of natural spline bases in each coordinate (total df = 4 x 4 = 16). The broken purple boundary is the Bayes decision boundary for this problem."   HTF

Multivariate Adaptive Regression Splines

Widely known by its acronym MARS. Apparently trademarked!

We consider the regression model, Y=f(X)+ϵ, where X∈(ℝ^p) and ϵ~NID(0,(σ^2)).  But extension of MARS are available for GLM.

Advantages

It is a fairly automatic, black-box procedure so perhaps it takes less skill.

Popular and well-known data mining method.  Supported by commerical software such as SAS and Salford Systems (http://www.salford-systems.com/products-mars.html) as well as being available in R and S-Plus.  In R we have the packages mda and earth which implement MARS.

Handles complex nonlinearity in moderately high dimensions.

Disadvantages

Considered too ad-hoc and not a proper statistical model. Too complex and ugly to understand.  Better to try to build a simpler model if that is possible.

If the goal is just prediction then MARS might be a good choice.

CRAN packages available: polspline and earth.

Jerome Friedman (http://www-stat.stanford.edu/~jhf/) invented MARS in a 1991 paper published in the Annals of Statistics. Originally he published fortran code to implement this method but it was removed when Salford Systems starting offering MARS software.  However, I have the original Fortran code.

Friedman,  J. H. (1991).  Multivariate Adaptive Regression Splines.  Annals of Statistics  19,  1–141.

The MARS Model

The model is built from simple linear spline functions, called reflected pairs,

(SubPlus[(x-t)])=(Piecewise[{{x-t, x>t}, {0, x≤t}}])

(SubPlus[(t-x)])=(Piecewise[{{t-x, x<t}, {0, x≥ t}}])

"Figure 9.9: The basis functions   (SubPlus[(x - t)])  (solid red) and (SubPlus[(t - x)])  (broken green) used by MARS."  HTF

Given n observations on p inputs, (Subscript[X, 1]),…,(Subscript[X, p]), denote the model matrix by 𝒳=(Subscript[((Subscript[x, i,j])), n×p]).

We start with the collection of basis functions

ℂ={(SubPlus[((Subscript[X, j])-t)]),(SubPlus[(t-(Subscript[X, j]))]): t∈{(Subscript[x, 1,j]),…,(Subscript[x, n,j])}, j=1,…,p}

The final regression model for, Y=f(X)+ϵ, can be written

f(X)=(Subscript[β, 0])+(∑_(m=1)^M)(Subscript[β, m])(Subscript[h, m])(X)

where (Subscript[h, m])(X) is comprised of functions in ℂ as well as products of two or functions in ℂ.  In the additive version of MARS attention product terms are excluded.  The additive version is the default in mda and earth but higher order degrees can be specified.  The coefficients (Subscript[β, 0]),(Subscript[β, 1]),…,(Subscript[β, m]) can be estimated by least squares.

MARS Model Fitting. First Stage.

In the first stage, we use forward stepwise regression to select the terms in the model.  Starting with (Subscript[h, 0])(X)=1, we find which reflected pair should enter next, to form {(Subscript[h, 0])(X), (Subscript[h, 1])(X),(Subscript[h, 2])(X)}, where (Subscript[h, 1]) and (Subscript[h, 2]) are of the form (Subscript[h, 1])(X)=(SubPlus[((Subscript[X, j])-t)]) and (Subscript[h, 2])(X)=(SubPlus[(t-(Subscript[X, j]))]) for some 1≤ j≤p for some t in the j-th column of X. At the next stage we consider all products of the form,

(Subscript[h, m])(X)(SubPlus[((Subscript[X, ℓ])-t)]) and (Subscript[h, m])(X)(SubPlus[(t-(Subscript[X, ℓ]))]) for all possible (Subscript[X, ℓ]),ℓ=1,…,p with t in the ℓ-th column of 𝒳.

In the additive case attention is restricted to, m=0, in each case.  Eventually the forward selection process stops after a large number of terms have been selected.

"Figure 9.10: Schematic of the MARS forward model building procedure. On the left are the basis functions currently in the model: initially, this is the constant function h(X) = 1. On the right are all candidate basis functions to be considered in building the model. These are pairs of piecewise linear basis functions as in Figure 9.9, with knots t at all unique observed values (Subscript[x, i,j]) of each predictor (Subscript[X, j]).. At each stage we consider all products of a candidate pair with a basis function in the model. The product that decreases the residual error the most is added into the current model. Above we illustrate the first three steps of the procedure, with the selected functions shown in red."  HTF

"Figure 9.11: The function h((Subscript[X, 1]) ,(Subscript[X, 2])) = (SubPlus[((Subscript[X, 1]) - (Subscript[x, 51]))])(SubPlus[((Subscript[x, 72]) - (Subscript[X, 2]))]) resulting from multiplication of two piecewise linear MARS basis functions."  HTF

MARS Model Fitting. Second Stage.

In the second stage of the fitting process, backward stepwise selection is used to find the best model, (Subscript[f, λ]), of size λ, where λ is the number of parameters estimated.  Because of its computational efficiency and stability GCV is used to choose among the best models of size λ,

GCV(λ)=(((∑_(i=1)^N)((((Subscript[y, i])-(Subscript[f, λ])((Subscript[x, i]))))^2))/(((1-M(λ)/n))^2))

where M(λ) is the effective number of parameters.  The effective number of parameters attempts to take into account the effect of knot selection.  In general M(λ)=r+c K, where r is the number of basis functions, K is the number of knots and c=2 in the additive case and c=3 in the non-additive case.

Spam Example

"Figure 9.12:  Spam data: test error misclassification rate jor the MARS procedure, as a junction oj the rank (number oj independent basis junctions) in the model." HTF

Mixture Classification Example

#Source: MARSMixture
#MARS: mixture training dataset
#Recall Bayes rate 20.76%
#Linear regression 27.3% (test data)
#Thin-plate splines 23.4%
library(mda)
library(earth)
d<-read.table("d:/math/2009/9850/data/TrainMixture.dat")
names(d)<-c("X","Y","Z")
n<-nrow(d)
#fit mars
ans<-mars(d[,1:2], d[,3], degree=2,newvar.penalty=1)
dnew<-read.table("d:/math/2009/9850/data/TestMixture.dat")
zHat<-predict(ans, newdata=dnew[,1:2])
zHat[zHat>=0.5]<-1
zHat[zHat<0.5]<-0
MisRate<-sum(zHat!=dnew[,3])/length(zHat)
MisRate

[1] 0.262

#fit earth - regular mars regression
ans<-earth(d[,1:2], d[,3], degree=2, newvar.penalty=1)
dnew<-read.table("d:/math/2009/9850/data/TestMixture.dat")
zHat<-predict(ans, newdata=dnew[,1:2])
zHat[zHat>=0.5]<-1
zHat[zHat<0.5]<-0
MisRate<-sum(zHat!=dnew[,3])/length(zHat)
MisRate

[1] 0.2485

#fit earth - logit mars regression
ans<-earth(d[,1:2], d[,3], degree=2, newvar.penalty=1, glm=list(family=binomial)) 
dnew<-read.table("d:/math/2009/9850/data/TestMixture.dat")
zHat<-predict(ans, newdata=dnew[,1:2])
zHat[zHat>=0.5]<-1
zHat[zHat<0.5]<-0
MisRate<-sum(zHat!=dnew[,3])/length(zHat)
MisRate

[1] 0.2685

Example with Sunspot Numbers

http://www.ngdc.noaa.gov/stp/SOLAR/ftpsunspotnumber.html

(Subscript[(OverHat[X]), t]) = 2.711 + .960(Subscript[X, t-1])+ .332(SubPlus[(47.0 - (Subscript[X, t-5]))])- .257(SubPlus[(59.1-(Subscript[X, t-9]))])- .003(SubPlus[(Subscript[X, t-1])((Subscript[X, t-2])- 26.0)])+ .017(SubPlus[(Subscript[X, t-1])(44.0 - (Subscript[X, t-3]))])- .032(SubPlus[(Subscript[X, t-1])(17.1 - (Subscript[X, t-4]))])+.004(SubPlus[(Subscript[X, t-1])(26.0 - (Subscript[X, t-2]))])(SubPlus[((Subscript[X, t-5]) - 41.0)])

"Table 4. The Mean Sum of Squares Error (Subscript[σ, ϵ]^2), Number of Model Parameters, and the Predictive Mean Sum of Squares Error (Subscript[σ, ϵ]^2)(i) for the ith Forward Step Prediction for the Period (1921-1955) of the AR, Bilinear, SETAR, and ASTAR Models of the Yearly Wolf Sunspot Numbers for the Period (1700-1920)."  Lewis (1991)

Lewis, P.A.  and Stevens,  J. G. (1991).  Nonlinear Modeling of Time Series Using Multivariate Adaptive Regression Splines,  JASA 86,  864-877.

Epidemiological time series

Many infectious diseases such as mumps and flu show chaotic seasonal behavior similar to the sunspot series such as the reported monthly cases of mumps, 1928 - 1972, New York City.

mumps={124,132,193,144,195,105,40,20,42,58,170,197,329,417,968,989,1264,1074,352,118,69,91,182,233,406,550,853,835,787,518,132,66,57,77,83,138,116,177,288,299,315,329,143,50,59,74,110,180,361,556,687,821,799,901,478,212,153,207,386,576,807,923,1604,1475,1123,662,253,116,70,107,147,155,266,249,406,452,547,472,213,136,114,137,276,341,688,935,1850,1938,1806,1030,472,151,123,133,195,255,357,386,543,645,668,659,276,142,93,82,122,225,315,360,777,1070,1115,1200,550,279,154,170,290,488,764,990,1695,1738,1295,1026,345,178,85,115,122,179,266,326,454,510,555,531,358,177,131,189,262,398,662,733,941,1229,1485,1406,959,436,233,221,391,696,930,899,1261,1245,847,626,348,168,132,142,201,399,569,651,1024,1272,1258,960,445,201,138,117,162,272,395,552,926,1062,1070,846,355,172,109,94,184,245,324,433,752,845,859,671,275,134,75,104,161,237,330,467,707,667,711,793,429,211,140,207,251,289,433,468,799,1138,1048,710,438,178,114,118,175,206,412,477,641,814,846,883,591,208,177,221,311,680,1105,1490,1956,1713,1291,1053,366,153,95,117,130,183,312,417,596,554,510,527,290,162,96,97,191,353,604,767,1116,1103,1330,1342,670,319,167,176,261,363,616,634,808,902,1003,833,475,243,167,152,182,279,520,615,745,838,787,827,570,254,182,185,268,530,813,785,1266,1495,1659,1532,1071,458,265,230,299,634,802,945,1220,1088,792,707,420,199,134,123,176,284,371,491,851,923,882,945,561,324,205,225,470,657,1066,1371,1674,1844,1819,1469,769,393,176,240,299,375,584,490,702,769,765,601,389,217,165,120,137,204,288,348,515,774,751,711,555,280,182,160,255,565,753,890,1183,1117,983,988,573,278,219,171,270,415,510,587,721,754,699,667,355,230,145,129,235,263,361,350,551,488,631,717,452,293,165,180,230,276,435,487,563,608,803,615,554,297,166,213,303,517,703,735,1037,1078,954,729,552,321,166,177,216,310,342,268,304,280,286,280,213,143,101,117,132,277,280,424,668,635,701,926,706,493,309,259,324,541,569,672,1020,848,754,765,437,255,166,130,165,166,235,261,390,424,500,440,294,167,126,124,162,182,274,349,420,607,572,496,415,240,158,124,135,100,164,148,298,506,567,639,447,283,158,223,212,333,330,281,354,527,463,429,300,147,104,90,77,122,146,149,248,300,223,294,235,110,88,68,88,132,149,157,219,221,264,298};
yr28=Transpose[{ConstantArray[1928,6],Range[7,12]}];
nyrs=1972-1927-1;
dts1=Transpose[{ConstantArray[#,12]&/@Range[1929,1972],Table[Range[1,12],{44}]},{3,1,2}];
dts=Join[{yr28},dts1];
dts0=Flatten[dts,1];
mumpsDts=Transpose[{dts0,mumps}];
DateListPlot[mumpsDts, FrameLabel->{"year","Number of Cases"},PlotLabel->"Reported monthly cases of mumps, 1928-1972, NYC"]

Data: Yorke, J.A. and London, W.P. (1973). “Recurrent Outbreaks of Measles, Chickenpox and Mumps”, American Journal of Epidemiology, Vol. 98, pp.469. There are several papers in the Chaos literature that model this series.

Simple Kriging and HM Derivation

We seek the minimum mean square error linear predictor for Y((Subscript[s, 0])), (Subscript[s, 0])∈(Subscript[𝔻, s]) based on the observed data Z((Subscript[s, i])),i=1,…,m .

CR (p. 137) provide a calculus based approach from first principles. The normal assumption is not needed. It may also be more simply derived using conditional mean and variance in a multivariate normal distribution.

Multivariate normal regression and variance

If we assume normality and for brievity write Y=Y((Subscript[s, 0])) and

Z=(({{Z((Subscript[s, 1]))}, {}, {Z((Subscript[s, m]))}}))

Then under the normal assumption,

(({{Y}, {Z}})) ~ Gau((({{(Subscript[μ, Y])}, {(Subscript[μ, Z])}})),(({{(Subscript[C, Y,Y]), (Subscript[C, Y,Z])}, {(Subscript[C, Z,Y]), (Subscript[C, Z,Z])}})))

the best linear predictor is given by the expected value,

𝔈{Y|Z}=(Subscript[μ, Y])+(Subscript[C, Y,Z])(Subscript[C, Z,Z]^(-1))(Z-(Subscript[μ, Z]))

and the prediction variance is

Var(Y|Z)=(Subscript[C, Y,Y])-(Subscript[C, Y,Z])(Subscript[C, Z,Z]^(-1))(Subscript[C, Z,Y])

Eqns. () and () give the general idea. But we need to give more explicit expressions for the terms involved. The derivation is completed using the HM formulation.

HM Formulation

Data Model

Conditional of (Subscript[σ, ε]^2) and taking i=1,…,m

Z((Subscript[s, i]))|Y((Subscript[s, i])) are IID Gaussian with mean (Subscript[μ, Z((Subscript[s, i]))])=Y((Subscript[s, i])) and variance (Subscript[σ, ε]^2).

Process Model

For i=1,…,m, and given p covariates x((Subscript[s, i]))=(((Subscript[x, 1])((Subscript[s, i])),…,(Subscript[x, p])((Subscript[s, p]))))' and regression coefficients, β,

Y((Subscript[s, i]))=(x((Subscript[s, i])))'β+u((Subscript[s, i]))

where u((Subscript[s, i])) is a spatial process with covariances defined by C(u((Subscript[s, i])),u((Subscript[s, j]))). So the covariance matrix of Y=((Y((Subscript[s, 1])),…,Y((Subscript[s, m]))))' is given by (Subscript[C, Y,Y])=((Subscript[C(u((Subscript[s, i])),u((Subscript[s, j]))), m×m]). The means are given by (Subscript[μ, Y])=(((x((Subscript[s, 1])))'β,…,(x((Subscript[s, m])))'β))'.

Expectation and Variance Conditional on Data

We want to find the conditional expectation and variance of Y((Subscript[s, 0])) given

Z=(({{Z((Subscript[s, 1]))}, {}, {Z((Subscript[s, m]))}}))

It is assumed that the parameters (Subscript[σ, ε]^2), μ((Subscript[s, i]))=(x((Subscript[s, i])))'β and C(u((Subscript[s, i]))), i=1,…,m are all known. Then the covariance matrix of Z is given by

Var(Z)=(Subscript[C, Z,Z])=(Subscript[C, Y,Y])+(Subscript[σ, ε]^2)(Subscript[ℐ, m])

and the expectation given the data is,

(Subscript[μ, Z])=X β

where

X=(({{(x((Subscript[s, 1])))'}, {⋮}, {x(((Subscript[s, m])))'}})).

Setting Y=Y((Subscript[s, 0])) with (Subscript[μ, Y])=(x((Subscript[s, 0])))'β, (Subscript[C, Y,Y])=C(Y((Subscript[s, 0])),Y((Subscript[s, 0])))≡(Subscript[σ, Y]^2), and

(Subscript[C, Y,Z])=(C(Y((Subscript[s, 0])),Z((Subscript[s, 1]))), …,C(Y((Subscript[s, 0])),Z((Subscript[s, m])))
=(C(Y((Subscript[s, 0])),Y((Subscript[s, 1]))), …,C(Y((Subscript[s, 0])),Y((Subscript[s, m])))

Hence we can write,

(({{Y}, {Z}})) ~ Gau((({{(Subscript[μ, Y])}, {(Subscript[μ, Z])}})),(({{(Subscript[C, Y,Y]), (Subscript[C, Y,Z])}, {(Subscript[C, Z,Y]), (Subscript[C, Z,Z])}})))

and the optimal prediction or Kriged estimate, (OverHat[Y])((Subscript[s, 0])), for Y((Subscript[s, 0])) is

(OverHat[Y])((Subscript[s, 0]))=(Subscript[μ, Y])+(Subscript[C, Y,Z])(Subscript[C, Z,Z]^(-1))(Z-(Subscript[μ, z]))
=(x((Subscript[s, 0])))'β+(C(Y((Subscript[s, 0])),Y((Subscript[s, 1]))), …,C(Y((Subscript[s, 0])),Y((Subscript[s, m]))) ((((Subscript[C, Y,Y])+(Subscript[σ, ε]^2)(Subscript[ℐ, m])))^(-1))(Z-X β)

and the variance is

Var((OverHat[Y])((Subscript[s, 0])))=(Subscript[σ, Y]^2)-C(Y((Subscript[s, 0])),Y((Subscript[s, m]))) ((((Subscript[C, Y,Y])+(Subscript[σ, ε]^2)(Subscript[ℐ, m])))^(-1))(C(Y((Subscript[s, 0])),Y((Subscript[s, m]))))'

Since (Subscript[C, Z,Z])=(Subscript[C, Y,Y])+(Subscript[σ, ε]^2)(Subscript[ℐ, m]) the predictor and its variance may also be written,

(OverHat[Y])((Subscript[s, 0]))=(x((Subscript[s, 0])))'β+(C(Y((Subscript[s, 0])),Y((Subscript[s, 1]))), …,C(Y((Subscript[s, 0])),Y((Subscript[s, m]))) (Subscript[C, Z,Z]^(-1))(Z-X β)

and

Var((OverHat[Y])((Subscript[s, 0])))=(Subscript[σ, Y]^2)-C(Y((Subscript[s, 0])),Y((Subscript[s, m])))  (Subscript[C, Z,Z]^(-1))(C(Y((Subscript[s, 0])),Y((Subscript[s, m]))))'

Decadal Temperature Differences over the Americas

The data a from a Climate Model that predicts monthly temperature. The resolution is 2.8°longitude by 2.8° latitude and data were obtained for a grid of 28×28 covering a portion of the Americas. Then the process was Y(s)=1990 to 1999 average - 1980 to 1989 average for each of the m=28×28=784 values. A heat map for the observed data is shown in the Figure below.

An exponential variogram was fitted to the detrended data, δ(s), obtained by using a polynomial or tensor spline in latitute and longitude, that is,

Y(s)=(x(s))'β+δ(s)

and shown by the blue curve. The red dots show the empirical variogram.

A 10×10 sub-grid was chosen and a Gaussian white noise random error was added to each value, for i=1,…,m; m=100,

Z((Subscript[s, i]))=Y((Subscript[s, i]))+ε((Subscript[s, i]))

The resulting grid is shown in the heatmap below.

The heat map below shows the predictions on the 28×28 grid using the simple Kriging predictor. It compares reasonably well to the original “truth”.

Kriging also provides estimates of the standard error of the predictions that are shown in the heat map below.

Non-Gaussian HM Kriging

Diggle, P.J., Tawn, J.A., and Moyeed, R.A. (1998). Model-based geostatistics. Applied Statistics 47, 299-350.

Diggle et al. (1998) develop an MCMC algorithm for Kriging using a GLM approach. The data Z is modelled using the exponential family. The process model Y is assumed Gaussian.

CW discuss some other non-Gaussian Kriging models.

Ordinary and Universal Kriging

These are the names given in the Classical Geostatistical Modelling approach.

The Universal Kriging Predictor Derivation

The Universal Kriging predictor, (OverBar[Y]), is derived as the unbiased linear predictor, λ'Z, with minimum variance and subject to the restriction (Subscript[λ, 1])+…+(Subscript[λ, m])=1. So it is like a weighted average, although some weights might be negative.

POINT FOR DISCUSSION. I find the discussion on CR that is repeated in his earlier book Cressie (1993) about homogenous functions vs non-homogenous functions obfusticating and unnecessary. So CR say eqn. () is a non-homogenous predictor because it is of the form λ'Z+k, where the k involves β. This is true but silly! Once the (OverHat[β]) is used, it becomes homogenous!! At first I thought there point was analgous to the question in time series ARMA modelling: whether or not to include a deterministic component such as a mean or linear trend.

The Universal Kriging Result

The derivations of these predictors is straightforward using calculus and the details are given CW (p. 148). The optimal universal predictor can be obtained from  eqn. () by replacing β by the GLS estimate,

(OverHat[β])=(((X'(Subscript[C, Z]^(-1))X))^(-1))X'(Subscript[C, Z]^(-1))Z.

The variance for (OverHat[Y])((Subscript[s, 0])) is slightly increased over than given in eqn. () due to estimation error for β since

Var((OverHat[β]))=(((X'(Subscript[C, Z]^(-1))X))^(-1)).

Remark. Computationally we would never use eqns. () and (). Matrix inverse is O((n^4)) complexity whereas solving linear equations is O((n^3)).  R uses the QR algorithm for LS computations.

The Ordinary Kriging Predictor

The Ordinary Kriging predictor is a special case of the Universal Kriging predictor assuming that 𝔈{Y((Subscript[s, i]))}=(Subscript[μ, Y]) does not depend any covariates or equivalently x((Subscript[s, i]))=1, i=1,…,m.

So, in other words, we estimate (Subscript[μ, Y]) using GLS estimator as in eqn. ().

Simple Kriging Predictor

𝔈{Y((Subscript[s, i]))}=(Subscript[μ, Y]) does not depend any covariates or equivalently x((Subscript[s, i]))=1, i=1,…,m. We estimate (Subscript[μ, Y]) using the sample mean.

BHM Kriging

The result is essentially unchanged when an improper uniform prior is put on β. It is still assumed that the covariance function is known or estimated from the data.

Total Column Ozone (TCO)

The original research for this was reported in 3 papers listed below in the References. LINK

Data

TOMS (Total Ozone Maping Spectrometer) data from October 5, 1998. Data resolution was 1° latitude and 1.25° longitude, resulting in 180×360/1.25=51840 pixels. This is (was) considered “massive datasets”. Regions that are white correspond to missing data. These are the observations Z, a vector of length 51840.

Results

EHM developed where Y latent variable of dimension 51840 based on the observations Z. The model parameters were estimated. The fixed rank Kriging method discussed by CW (§4.14, eqns. 4.80 and 4.81) provides estimates of Y((Subscript[s, 0])) and its covariance matrix with complexity only O(n) where as the usual GLS based methods require O((n^3)).

Posterior mean of Y in panel (a) and its standard deviation in panel (b).

Large and Massive Datasets

My current PC has only 16GB RAM and it easily handles matricies as large as 5184. Complexity of matrix inverse is O((n^4)) in the general case. Presumably Mathematica does not detect that the matrix is in fact a Toeplitz and special methods that have O((n^3)) complexity are available (McLeod, Yu and Kroughly 2005).

With[{n=5184},
Timing[Inverse[ToeplitzMatrix[(1/Range[n])^0.25]];]]

But it is not big enough to invert the above matrix when n=51840. Sometime I will try with my new linux laptop that has 32GB.

References

[JCH] Johannsesson, Cressie, and Huang (2007).  Dynamic multi-resolution spatial models. Environmental and Ecological Statistics, 14/5 5-25.

[HCG] Huang, Cressie and Gabrosek (2002). Fast resolution-consistent spatial prediction of global processes from satellite data. Journal of Computational Statistics and Graphics, 11 63-88.

[THC] Tzeng, Huang, and Cressie (2005). A fast optimal spatial-prediction method for massive datasets. Journal of the American Statistical Association 100 1343-1357.