autoFRK: Automatic Fixed Rank Kriging
In autoFRK: Automatic Fixed Rank Kriging

Description Usage Arguments Details Value Author(s) References See Also Examples

This function performs resolution adaptive fixed rank kriging based on spatial data observed at one or multiple time points via the following spatial random-effects model:

z[t]=μ +G*w[t]+η[t]+e[t], w[t]~N(0,M); e[t] ~ N(0, s*D); t=1,...,T,

where z[t] is an n-vector of (partially) observed data at n locations, μ is an n-vector of deterministic mean values, D is a given n by n matrix, G is a given n by K matrix, η[t] is an n-vector of random variables corresponding to a spatial stationary process, and w[t] is a K-vector of unobservable random weights. Parameters are estimated by maximum likelihood in a closed-form expression. The matrix G corresponding to basis functions is given by an ordered class of thin-plate spline functions, with the number of basis functions selected by Akaike's information criterion.

autoFRK(
  Data,
  loc,
  mu = 0,
  D = diag.spam(NROW(Data)),
  G = NULL,
  finescale = FALSE,
  maxit = 50,
  tolerance = 0.1^6,
  maxK = NULL,
  Kseq = NULL,
  method = c("fast", "EM"),
  n.neighbor = 3,
  maxknot = 5000
)

`Data`	n by T data matrix (NA allowed) with z[t] as the t-th column.
`loc`	n by d matrix of coordinates corresponding to n locations.
`mu`	n-vector or scalar for μ; Default is 0.
`D`	n by n matrix (preferably sparse) for the covariance matrix of the measurement errors up to a constant scale. Default is an identity matrix.
`G`	n by K matrix of basis function values with each column being a basis function taken values at `loc`. Default is NULL, which is automatic determined.
`finescale`	logical; if `TRUE` then a (approximate) stationary finer scale process η[t] will be included based on `LatticeKrig` pacakge. In such a case, only the diagonals of D would be taken into account. Default is `FALSE`.
`maxit`	maximum number of iterations. Default is 50.
`tolerance`	precision tolerance for convergence check. Default is 0.1^6.
`maxK`	maximum number of basis functions considered. Default is 10 \cdot √{n} for n>100 or n for n<=100.
`Kseq`	user-specified vector of numbers of basis functions considered. Default is `NULL`, which is determined from `maxK`.
`method`	"fast" or " EM"; if "fast" then the missing data are filled in using k-nearest-neighbor imputation; if "EM" then the missing data are taken care by the EM algorithm. Default is "fast".
`n.neighbor`	number of neighbors to be used in the "fast" imputation method. Default is 3.
`maxknot`	maximum number of knots to be used in generating basis functions. Default is 5000.

The function computes the ML estimate of M using the closed-form expression in Tzeng and Huang (2018). If the user would like to specify a D other than an identity matrix for a large n, it is better to provided via spam function in spam package.

an object of class FRK is returned, which is a list of the following components:

`M`	ML estimate of M.
`s`	estimate for the scale parameter of measurement errors.
`negloglik`	negative log-likelihood.
`w`	K by T matrix with w[t] as the t-th column.
`V`	K by K matrix of the prediction error covariance matrix of w[t].
`G`	user specified basis function matrix or an automatically generated `mrts` object.
`LKobj`	a list from calling `LKrig.MLE` in `LatticeKrig` package if `useLK=TRUE`; otherwise `NULL`. See that package for details.

ShengLi Tzeng and Hsin-Cheng Huang.

Tzeng, S., & Huang, H.-C. (2018). Resolution Adaptive Fixed Rank Kriging, Technometrics, https://doi.org/10.1080/00401706.2017.1345701.

predict.FRK

#### generating data from two eigenfunctions
originalPar <- par(no.readonly = TRUE)
set.seed(0)
n <- 150
s <- 5
grid1 <- grid2 <- seq(0, 1, l = 30)
grids <- expand.grid(grid1, grid2)
fn <- matrix(0, 900, 2)
fn[, 1] <- cos(sqrt((grids[, 1] - 0)^2 + (grids[, 2] - 1)^2) * pi)
fn[, 2] <- cos(sqrt((grids[, 1] - 0.75)^2 + (grids[, 2] - 0.25)^2) * 2 * pi)

#### single realization simulation example
w <- c(rnorm(1, sd = 5), rnorm(1, sd = 3))
y <- fn %*% w
obs <- sample(900, n)
z <- y[obs] + rnorm(n) * sqrt(s)
X <- grids[obs, ]

#### method1: automatic selection and prediction
one.imat <- autoFRK(Data = z, loc = X, maxK = 15)
yhat <- predict(one.imat, newloc = grids)

#### method2: user-specified basis functions
G <- mrts(X, 15)
Gpred <- predict(G, newx = grids)
one.usr <- autoFRK(Data = z, loc = X, G = G)
yhat2 <- predict(one.usr, newloc = grids, basis = Gpred)

require(fields)
par(mfrow = c(2, 2))
image.plot(matrix(y, 30, 30), main = "True")
points(X, cex = 0.5, col = "grey")
image.plot(matrix(yhat$pred.value, 30, 30), main = "Predicted")
points(X, cex = 0.5, col = "grey")
image.plot(matrix(yhat2$pred.value, 30, 30), main = "Predicted (method 2)")
points(X, cex = 0.5, col = "grey")
plot(yhat$pred.value, yhat2$pred.value, mgp = c(2, 0.5, 0))
par(originalPar)
#### end of single realization simulation example

#### independent multi-realization simulation example
set.seed(0)
wt <- matrix(0, 2, 20)
for (tt in 1:20) wt[, tt] <- c(rnorm(1, sd = 5), rnorm(1, sd = 3))
yt <- fn %*% wt
obs <- sample(900, n)
zt <- yt[obs, ] + matrix(rnorm(n * 20), n, 20) * sqrt(s)
X <- grids[obs, ]
multi.imat <- autoFRK(Data = zt, loc = X, maxK = 15)
Gpred <- predict(multi.imat$G, newx = grids)

G <- multi.imat$G
Mhat <- multi.imat$M
dec <- eigen(G %*% Mhat %*% t(G))
fhat <- Gpred %*% Mhat %*% t(G) %*% dec$vector[, 1:2]
par(mfrow = c(2, 2))
image.plot(matrix(fn[, 1], 30, 30), main = "True Eigenfn 1")
image.plot(matrix(fn[, 2], 30, 30), main = "True Eigenfn 2")
image.plot(matrix(fhat[, 1], 30, 30), main = "Estimated Eigenfn 1")
image.plot(matrix(fhat[, 2], 30, 30), main = "Estimated Eigenfn 2")
par(originalPar)
#### end of independent multi-realization simulation example