Title: | Tools for Spatial Data Analysis |
---|---|
Description: | Tools for spatial data analysis. Emphasis on kriging. Provides functions for prediction and simulation. Intended to be relatively straightforward, fast, and flexible. |
Authors: | Joshua French <[email protected]> |
Maintainer: | Joshua French <[email protected]> |
License: | GPL (>=2) |
Version: | 1.0.5 |
Built: | 2024-11-06 04:41:46 UTC |
Source: | https://github.com/jfrench/spatialtools |
coincident
takes the coordinate matrices
coords1
and coords2
and returns the indices
of the coincident coordinates of the two matrices.
coincident(coords1, coords2)
coincident(coords1, coords2)
coords1 |
An |
coords2 |
An |
This function calls a compiled C++
program created
using the Rcpp package. This (may) result in a significant
speedup over pure R code since the search algorithm
involves loops. We assume that there are no duplicate
coordinates in coords1
and coords2
,
respectively. In other words, each row of coords1
is
unique and each row of coords2
is unique. There is
at most 1 row of coords1
that will match with a row
in coords2
.
Returns a matrix with the indices of the coincident
locations. Specifically, an matrix will be
returned, with
being the number of coordinates in
coords1
coinciding with coordinates in
coords2
. If row i
of the matrix is c(2, 5),
then the i
th set of coincident locations is between
the 2nd row of coords1
and the 5th row of
coords2
. If there are no coincident locations, then
a matrix of size is returned.
Joshua French
#Generate two sets of coordinates loc1 <- as.matrix(expand.grid(seq(0, 1, len = 25), seq(0, 1, len = 25))) loc2 <- as.matrix(expand.grid(seq(0, 1, len = 101), seq(0, 1, len = 101))) coincident(loc1, loc2)
#Generate two sets of coordinates loc1 <- as.matrix(expand.grid(seq(0, 1, len = 25), seq(0, 1, len = 25))) loc2 <- as.matrix(expand.grid(seq(0, 1, len = 101), seq(0, 1, len = 101))) coincident(loc1, loc2)
Calculates spatial covariance matrix of the observed responses, and
possibly, the responses to be predicted.
If pcoords
is not provided, then only V
,
the covariance matrix of the observed responses will be returned.
If pcoords
is provided, then Vp
and Vop
(the covariance matrix for predicted responses and between observed and
predicted responses, respectively) will also be returned.
cov.sp(coords, sp.type = "exponential", sp.par = stop("specify sp.par argument"), error.var = 0, smoothness = 0.5, finescale.var = 0, pcoords = NULL, D = NULL, Dp = NULL, Dop = NULL)
cov.sp(coords, sp.type = "exponential", sp.par = stop("specify sp.par argument"), error.var = 0, smoothness = 0.5, finescale.var = 0, pcoords = NULL, D = NULL, Dp = NULL, Dop = NULL)
coords |
A numeric matrix of size |
sp.type |
A character vector specifying the spatial covariance type. Valid types are currently exponential, gaussian, matern, matern2, and spherical. |
sp.par |
A vector of length 2 specifying the scale and strength of dependence of the covariance function. The first element is the variance of the underlying spatial process (also known as the hidden or latent spatial process). This value is also called the partial sill. The second element is the strength of dependence between responses. |
error.var |
A non-negative number indicating the variance of the error term. |
smoothness |
A positive number indicating the variance of the error term. |
finescale.var |
A non-negative positive number indicating the finescale variability. The is also called the microscale variance |
pcoords |
A numeric matrix of size |
D |
The Euclidean distance matrix for the |
Dp |
The Eucliean distance matrix for the |
Dop |
The Euclidean intersite distance matrix between the locations in |
The spatial covariance functions are parameterized in a manner
consistent with the cov.spatial
function in the
geoR
package. The matern2
covariance function is an alternative covariance function suggested by Handcock and Wallis (1994).
The benefit of this parameterization is that the range parameter is that it allows the effective range to be less dependent on the smoothness parameter.
The D
, Dp
, and Dop
arguments are supplied to decrease the number of necessary computations needed when performing repetitive analysis or simulations. It is probably in the user's interest to not supply these arguments unless the duration of analysis is an important consideration. Note that these arguments override the information given in coords
and pcoords
, i.e., if dist1(coords) != D, then D is used in subsequent calculations, etc. This could create problems.
Returns a list with the following elements:
V |
The covariance matrix for the observed responses. Will be of size |
Vp |
The covariance matrix for the predicted responses. Only returned if |
Vp |
The covariance matrix between the observed responses and the predicted responses. Only returned if |
Joshua French
M.S. Handcock, J.R. Wallis. An approach to statistical spatial-temporal modeling of meteorological fields (with discussion). Journal of the American Statistical Association, 89 (1994), pp. 368–390.
simple.cov.sp
coords <- matrix(rnorm(30), ncol = 3) cov.sp(coords = coords, sp.type = "exponential", sp.par = c(2, 1), error.var = 1)
coords <- matrix(rnorm(30), ncol = 3) cov.sp(coords = coords, sp.type = "exponential", sp.par = c(2, 1), error.var = 1)
Calculates spatial covariance matrix of the observed responses, and
possibly, the responses to be predicted.
If pcoords
is not provided, then only V
,
the covariance matrix of the observed responses will be returned.
If pcoords
is provided, then Vp
and Vop
(the covariance matrix for predicted responses and between observed and
predicted responses, respectively) will also be returned.
cov.st(coords, time, sp.type = "exponential", sp.par = stop("specify sp.par argument"), error.var = 0, smoothness = 0.5, finescale.var = 0, t.type = "ar1", t.par = .5, pcoords = NULL, ptime = NULL, D = NULL, Dp = NULL, Dop = NULL, T = NULL, Tp = NULL, Top = NULL)
cov.st(coords, time, sp.type = "exponential", sp.par = stop("specify sp.par argument"), error.var = 0, smoothness = 0.5, finescale.var = 0, t.type = "ar1", t.par = .5, pcoords = NULL, ptime = NULL, D = NULL, Dp = NULL, Dop = NULL, T = NULL, Tp = NULL, Top = NULL)
coords |
A numeric matrix of size |
time |
A numeric matrix of size |
sp.type |
A character vector specifying the spatial covariance type. Valid types are currently exponential, gaussian, matern, and spherical. |
sp.par |
A vector of length 2 specifying the scale and dependence of the covariance function. The first element refers to the variance of the hidden process (sometimes this is called the partial sill) while the second elements determines the strength of dependence between locations. |
error.var |
A non-negative number indicating the variance of the error term. |
smoothness |
A positive number indicating the variance of the error term. |
finescale.var |
A non-negative positive number indicating the finescale variability. The is also called the microscale variance |
t.type |
A character vector indicating the temporal dependance structure. Currently, only "ar1" is implemented. |
t.par |
A numeric vector of length 1 indicating the strength of temporal dependence. |
pcoords |
A numeric matrix of size |
ptime |
A numeric matrix of size |
D |
The Euclidean distance matrix for the coords matrix. Must be of size |
Dp |
The Eucliean distance matrix for the pcoords matrix. Must be of size |
Dop |
The Eucliean intersite distance matrix between the locations in coords and the locations in pcoords. Must be of size |
T |
The Euclidean distance matrix for the time matrix. Must be of size |
Tp |
The Eucliean distance matrix for the ptime matrix. Must be of size |
Top |
The Eucliean intertime distance matrix between the times in |
At this point, this function only implements a separable spatio-temporal covariance funcation. If is the distance between
two sites, and
is the temporal lag between the times when the associated responses were observed, then the covariance function
where
is a spatial covariance function corresponding to the
exponential, matern, gaussian, or spherical and
is the temporal covariance function corresponding to an ar1 process with
.
The D
, Dp
, Dop
, T
, Tp
, Top
arguments are supplied to decrease the number of necessary computations needed when performing repetitive analysis or simulations. It is probably in the user's interest to not supply these arguments unless the duration of analysis is an important consideration. Note that these arguments override the information given in coords
, pcoords
, time
, and prime
, i.e., if dist1(coords) != D, then D is used in subsequent calculations, etc. This could create problems.
Returns a list with the following elements:
V |
The covariance matrix for the observed responses. |
Vp |
The covariance matrix for the predicted responses. Only returned if |
Vp |
The covariance matrix between the observed responses and the predicted responses. Only returned if |
Joshua French
simple.cov.sp
coords <- matrix(rnorm(30), ncol = 3) pcoords <- matrix(rnorm(90), ncol = 3) time <- 1:10 ptime <- 1:30 cov.st(coords = coords, time = time, sp.type = "exponential", sp.par = c(2, 1), error.var = 1, t.type = "ar1", t.par = .5, pcoords = pcoords, ptime = ptime)
coords <- matrix(rnorm(30), ncol = 3) pcoords <- matrix(rnorm(90), ncol = 3) time <- 1:10 ptime <- 1:30 cov.st(coords = coords, time = time, sp.type = "exponential", sp.par = c(2, 1), error.var = 1, t.type = "ar1", t.par = .5, pcoords = pcoords, ptime = ptime)
Calculates a decomposition of the provided covariance matrix, V
, using the chosen method.
decomp.cov(V, method = "eigen")
decomp.cov(V, method = "eigen")
V |
A (symmetric, positive-definite) covariance matrix. |
method |
A character vector specifying the method used to decompose |
The matrix V
is assumed to be symmetric and positive definite. Symmetry is checked, but the positive definiteness of the matrix is not. Returns a decomposition matrix U
such that V
= U
%*% t(U)
.
Returns a decomposition matrix U
such that V
= U
%*% t(U)
.
Joshua French
cov.sp
data(toydata) U <- decomp.cov(toydata$V, method = "chol") #range(toydata$V - U %*% t(U))
data(toydata) U <- decomp.cov(toydata$V, method = "chol") #range(toydata$V - U %*% t(U))
dist1
takes a matrix of coordinates and returns the
Euclidean distance matrix of the coordinates. It does this
using a compiled C
program, so it is faster than the
builtin R dist
function.
dist1(coords)
dist1(coords)
coords |
An |
An matrix of Euclidean distances.
Joshua French
x <- matrix(rnorm(30), ncol = 3) dist1(x)
x <- matrix(rnorm(30), ncol = 3) dist1(x)
dist2
takes the matrices of coordinates
coords1
and coords2
and returns the
inter-Euclidean distances between coordinates.
dist2(coords1, coords2)
dist2(coords1, coords2)
coords1 |
An |
coords2 |
An |
An matrix of Euclidean distances.
Joshua French
dist, dist1
x1 <- matrix(rnorm(30), ncol = 3) x2 <- matrix(rnorm(60), ncol = 3) dist2(x1, x2)
x1 <- matrix(rnorm(30), ncol = 3) x2 <- matrix(rnorm(60), ncol = 3) dist2(x1, x2)
Takes contours of contourLines
function and extracts the associated coordinates.
get.contours(x)
get.contours(x)
x |
A list returned by the |
Returns a 2-column matrix containing the coordinates making up the contours in contours.list.
Joshua French
contourLines, contour
data(volcano) x <- 10*1:nrow(volcano) y <- 10*1:ncol(volcano) cL <- contourLines(x, y, volcano) out <- get.contours(cL)
data(volcano) x <- 10*1:nrow(volcano) y <- 10*1:ncol(volcano) cL <- contourLines(x, y, volcano) out <- get.contours(cL)
Performs Ordinary Kriging using y
, the matrix of observed responses,
V
, the (positive definite) covariance matrix of the
observed responses, Vp
, the
covariance matrix of the responses to be predicted, and
Vop
,
the matrix of covariances between the observed
responses and the responses to be predicted.
krige.ok(y, V, Vp, Vop, nsim = 0, Ve.diag = NULL, method = "eigen")
krige.ok(y, V, Vp, Vop, nsim = 0, Ve.diag = NULL, method = "eigen")
y |
The vector of observed responses.
Should be a matrix of size |
V |
The covariance matrix of the observed responses.
The size is |
Vp |
The covariance matrix of the responses to be predicted.
The size is |
Vop |
The cross-covariance between the observed responses
and the responses to be predicted. The size is
|
nsim |
The number of simulated data sets to sample from the conditional predictive distribution. |
Ve.diag |
A vector of length |
method |
The method for decomposing |
It is assumed that there are observed data values
and that we wish to make predictions at
locations.
If doing conditional simulation, the Cholesky decomposition should not work when there are coincident locations between the observed data locations and the predicted data locations. Both the Eigen and Singular Value Decompositions should work.
If user specifies nsim
to be a positive integer, then nsim
conditional realizations of the predictive distribution will be generated. If this is less than 1, then no conditional simulation is done. If nsim
is a positive integer, then Ve.diag
must also be supplied. Ve.diag
is should be a vector of length specifying the measurement error variances of the observed data. This information is only used for conditional simulation, so this argument is only needed when
nsim
> 0. When conditional simulation is desired, then the argument method
can be to specify the method used to decompose V
. Options are "eigen", "chol", or "svd" (Eigen decomposition, Cholesky decomposition, or Singular value decomposition, respectively). This information is only used for conditional simulation, so this argument is only applicable when nsim
> 0.
The function returns a list containing the following objects:
pred |
A vector of length |
mspe |
A vector of length |
coeff |
A vector of length |
vcov.coeff |
A |
simulations |
An |
If nsim
> 0, this list has class "krigeConditionalSample".
Joshua French
Statistical Methods for Spatial Data Analysis, Schabenberger and Gotway (2003). See p. 226-228.
# create observed and predicted coordinates ocoords <- matrix(runif(100), ncol = 2) pcoords <- matrix(runif(200), ncol = 2) # include some observed locations in the predicted coordinates acoords <- rbind(ocoords, pcoords) # create covariance matrix C3 <- cov.sp(coords = ocoords, sp.type = "matern", sp.par = c(2, 1), smoothness = 1, finescale = 0, error = 0.5, pcoords = acoords) # generate data with error y <- rmvnorm(nsim = 1, mu = rep(2, 50), V = C3$V) + rnorm(50, sd = sqrt(.5)) # use universal kriging to make predictions. Do not do conditional simulation krige.obj <- krige.ok(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, nsim = 0) #Do conditional simulation krige.obj2 <- krige.ok(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, nsim = 100, Ve.diag = rep(.5, 50), method = "eigen")
# create observed and predicted coordinates ocoords <- matrix(runif(100), ncol = 2) pcoords <- matrix(runif(200), ncol = 2) # include some observed locations in the predicted coordinates acoords <- rbind(ocoords, pcoords) # create covariance matrix C3 <- cov.sp(coords = ocoords, sp.type = "matern", sp.par = c(2, 1), smoothness = 1, finescale = 0, error = 0.5, pcoords = acoords) # generate data with error y <- rmvnorm(nsim = 1, mu = rep(2, 50), V = C3$V) + rnorm(50, sd = sqrt(.5)) # use universal kriging to make predictions. Do not do conditional simulation krige.obj <- krige.ok(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, nsim = 0) #Do conditional simulation krige.obj2 <- krige.ok(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, nsim = 100, Ve.diag = rep(.5, 50), method = "eigen")
Performs simple kriging using y
, a vector of length ,
V
, the (positive definite) covariance matrix of the
observed responses, Vp
, the
covariance matrix of the responses to be predicted,
Vop
,
the matrix of covariances between the observed
responses and the responses to be predicted, and
m
, a numeric vector
of length 1 identifying the value of the mean
for each response.
krige.sk(y, V, Vp, Vop, m = 0, nsim = 0, Ve.diag = NULL, method = "eigen")
krige.sk(y, V, Vp, Vop, m = 0, nsim = 0, Ve.diag = NULL, method = "eigen")
y |
The vector of observed responses.
Should be a matrix of size |
V |
The covariance matrix of the observed responses.
The size is |
Vp |
The covariance matrix of the responses to be predicted.
The size is |
Vop |
The cross-covariance between the observed responses
and the responses to be predicted. The size is
|
m |
A numeric vector of length 1 giving the mean of each response. |
nsim |
The number of simulated data sets to sample from the conditional predictive distribution. |
Ve.diag |
A vector of length |
method |
The method for decomposing |
It is assumed that there are observed data values
and that we wish to make predictions at
locations.
The mean is subtracted from each value of
y
before determining the kriging weights,
and then the mean is added onto the predicted response.
If doing conditional simulation, the Cholesky decomposition should not work when there are coincident locations between the observed data locations and the predicted data locations. Both the Eigen and Singular Value Decompositions should work.
If user specifies nsim
to be a positive integer, then nsim
conditional realizations of the predictive distribution will be generated. If this is less than 1, then no conditional simulation is done. If nsim
is a positive integer, then Ve.diag
must also be supplied. Ve.diag
is should be a vector of length specifying the measurement error variances of the observed data. This information is only used for conditional simulation, so this argument is only needed when
nsim
> 0. When conditional simulation is desired, then the argument method
can be to specify the method used to decompose V
. Options are "eigen", "chol", or "svd" (Eigen decomposition, Cholesky decomposition, or Singular value decomposition, respectively). This information is only used for conditional simulation, so this argument is only applicable when nsim
> 0.
The function returns a list containing the following objects:
pred |
A vector of length |
mspe |
A vector of length |
simulations |
An |
mean |
The mean value (m) originally provided to the function |
.
If nsim
> 0, this list has class "krigeConditionalSample".
Joshua French
Statistical Methods for Spatial Data Analysis, Schabenberger and Gotway (2003). See p. 226-228.
data(toydata) y <- as.vector(toydata$y) V <- toydata$V Vp <- toydata$Vp Vop <- toydata$Vop krige.sk(y, V, Vp, Vop, m = 2)
data(toydata) y <- as.vector(toydata$y) V <- toydata$V Vp <- toydata$Vp Vop <- toydata$Vop krige.sk(y, V, Vp, Vop, m = 2)
Performs universal kriging using X
, the
design matrix for the regression coefficients of the observed
data,
y
, the matrix of observed responses,
V
, the (positive definite) covariance matrix of the
observed responses, Xp
, the design matrix
of the responses to be predicted,
Vp
, the
covariance matrix of the responses to be predicted, and
Vop
,
the matrix of covariances between the observed
responses and the responses to be predicted. If user specifies
nsim
to be a positive integer, then nsim
conditional realizations of the predictive distribution will be generated.
krige.uk(y, V, Vp, Vop, X, Xp, nsim = 0, Ve.diag = NULL, method = "eigen")
krige.uk(y, V, Vp, Vop, X, Xp, nsim = 0, Ve.diag = NULL, method = "eigen")
y |
The vector of observed responses. Should be a matrix of size |
V |
The covariance matrix of the observed responses.
The size is |
Vp |
The covariance matrix of the responses to be predicted. The size is |
Vop |
The cross-covariance between the observed responses and the responses to be predicted. The size is |
X |
The design matrix of the observed data. The size is |
Xp |
The design matrix of the responses to be predicted.
The size is |
.
nsim |
The number of simulated data sets to sample from the conditional predictive distribution. |
Ve.diag |
A vector of length |
method |
The method for decomposing |
It is assumed that there are observed data values
and that we wish to make predictions at
locations. We assume
that there are
regression coefficients (including the intercept).
Both
X
and Xp
should contain a column of 1's if an intercept
is desired.
If doing conditional simulation, the Cholesky decomposition should not work when there are coincident locations between the observed data locations and the predicted data locations. Both the Eigen and Singular Value Decompositions should work.
If user specifies nsim
to be a positive integer, then nsim
conditional realizations of the predictive distribution will be generated. If this is less than 1, then no conditional simulation is done. If nsim
is a positive integer, then Ve.diag
must also be supplied. Ve.diag
is should be a vector of length specifying the measurement error variances of the observed data. This information is only used for conditional simulation, so this argument is only needed when
nsim
> 0. When conditional simulation is desired, then the argument method
can be to specify the method used to decompose V
. Options are "eigen", "chol", or "svd" (Eigen decomposition, Cholesky decomposition, or Singular value decomposition, respectively). This information is only used for conditional simulation, so this argument is only applicable when nsim
> 0.
The function returns a list containing the following objects:
pred |
A vector of length |
mspe |
A vector of length |
coeff |
A vector of length |
vcov.coeff |
A |
sim |
An |
If nsim
> 0, this list has class "krigeConditionalSample".
Joshua French
Statistical Methods for Spatial Data Analysis, Schabenberger and Gotway (2003). See p. 241-243.
# create observed and predicted coordinates ocoords <- matrix(runif(100), ncol = 2) pcoords <- matrix(runif(200), ncol = 2) # include some observed locations in the predicted coordinates acoords <- rbind(ocoords, pcoords) # create design matrices X <- as.matrix(cbind(1, ocoords)) Xa <- as.matrix(cbind(1, acoords)) # create covariance matrix C3 <- cov.sp(coords = ocoords, sp.type = "matern", sp.par = c(2, 1), smoothness = 1, finescale = 0, error = 0.5, pcoords = acoords) # set values of regression coefficients coeff <- matrix(c(1, 2, 3), nrow = 1) # generate data with error y <- rmvnorm(nsim = 1, mu = tcrossprod(X, coeff), V = C3$V) + rnorm(50, sd = sqrt(.5)) # use universal kriging to make predictions. Do not do # conditional simulation krige.obj <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, X = X, Xp = Xa, nsim = 0) #Do kriging with conditional simulation krige.obj2 <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, X = X, Xp = Xa, nsim = 100, Ve.diag = rep(.5, 50), method = "eigen")
# create observed and predicted coordinates ocoords <- matrix(runif(100), ncol = 2) pcoords <- matrix(runif(200), ncol = 2) # include some observed locations in the predicted coordinates acoords <- rbind(ocoords, pcoords) # create design matrices X <- as.matrix(cbind(1, ocoords)) Xa <- as.matrix(cbind(1, acoords)) # create covariance matrix C3 <- cov.sp(coords = ocoords, sp.type = "matern", sp.par = c(2, 1), smoothness = 1, finescale = 0, error = 0.5, pcoords = acoords) # set values of regression coefficients coeff <- matrix(c(1, 2, 3), nrow = 1) # generate data with error y <- rmvnorm(nsim = 1, mu = tcrossprod(X, coeff), V = C3$V) + rnorm(50, sd = sqrt(.5)) # use universal kriging to make predictions. Do not do # conditional simulation krige.obj <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, X = X, Xp = Xa, nsim = 0) #Do kriging with conditional simulation krige.obj2 <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, X = X, Xp = Xa, nsim = 100, Ve.diag = rep(.5, 50), method = "eigen")
Estimates covariance parameters of spatial covariance functions using maximum likelihood or restricted maximum likelihood. See cov.sp
for more details of covariance functions to be estimated.
maxlik.cov.sp(X, y, coords, sp.type = "exponential", range.par = stop("specify range.par argument"), error.ratio = stop("specify error.ratio argument"), smoothness = 0.5, D = NULL, reml = TRUE, lower = NULL, upper = NULL, control = list(trace = TRUE), optimizer="nlminb")
maxlik.cov.sp(X, y, coords, sp.type = "exponential", range.par = stop("specify range.par argument"), error.ratio = stop("specify error.ratio argument"), smoothness = 0.5, D = NULL, reml = TRUE, lower = NULL, upper = NULL, control = list(trace = TRUE), optimizer="nlminb")
X |
A numeric matrix of size |
y |
A vector of length |
coords |
A numeric matrix of size |
sp.type |
A character vector specifying the spatial covariance type. Valid types are currently exponential, gaussian, matern, and spherical. |
range.par |
An initial guess for the spatial dependence parameter. |
error.ratio |
A value non-negative value indicating the ratio |
smoothness |
A positive number indicating the smoothness of the matern covariance function, if applicable. |
D |
The Euclidean distance matrix for the coords matrix. Must be of size |
reml |
A boolean value indicating whether restricted maximum likelihood estimation should be used. Defaults to TRUE. |
lower |
A vector giving lower bounds for the covariance parameters |
upper |
A vector giving upper bounds for the covariance parameters |
control |
A list giving tuning parameters for the |
optimizer |
A vector describing the optimization function to use for the optimization. Currently, only |
When doing the numerical optimizaiton, the covariance function is reparameterized slightly to speedup computation.
Specifically, the variance parameter for the process of interest,sp.par[1]
, is profiled out,
and the error.var
parameter is parameterized as sp.par[1] * error.ratio
, where error.ratio = error.var/sp.par[1]
.
Returns a list with the following elements:
sp.type |
The covariance form used. |
sp.par |
A vector containing the estimated variance of the hidden process and the spatial dependence. |
error.var |
The estimated error variance. |
smoothness |
The smoothness of the matern covariance function. |
par |
The final values of the optimization parameters. Note that these will not necessarily match |
convergence |
Convergence message from |
message |
Message from |
iterations |
Number of iterations for optimization to converge. |
evaluations |
Evaluations from |
Joshua French
cov.st
#generate 20 random (x, y) coordinates coords <- matrix(rnorm(20), ncol = 2) #create design matrix X <- cbind(1, coords) #create mean for observed data to be generated mu <- X %*% c(1, 2, 3) #generate covariance matrix V <- exp(-dist1(coords)) #generate observe data y <- rmvnorm(mu = mu, V = V) #find maximum likelihood estimates of covariance parameters maxlik.cov.sp(X = X, y = y, coords = coords, sp.type = "exponential", range.par = 1, error.ratio = 0, reml = TRUE)
#generate 20 random (x, y) coordinates coords <- matrix(rnorm(20), ncol = 2) #create design matrix X <- cbind(1, coords) #create mean for observed data to be generated mu <- X %*% c(1, 2, 3) #generate covariance matrix V <- exp(-dist1(coords)) #generate observe data y <- rmvnorm(mu = mu, V = V) #find maximum likelihood estimates of covariance parameters maxlik.cov.sp(X = X, y = y, coords = coords, sp.type = "exponential", range.par = 1, error.ratio = 0, reml = TRUE)
Estimates covariance parameters of spatio-temporal covariance functions using maximum likelihood or restricted maximum likelihood. See cov.st
for more details of covariance functions to be estimated. The covariance function is reparameterized slightly to speedup computation. Specifically, the variance parameter for the hidden process, sp.par[1], is profiled out and the error.var parameter is parameterized as sp.par[1] * error.ratio.
maxlik.cov.st(X, y, coords, time, sp.type = "exponential", range.par = stop("specify range.par argument"), error.ratio = stop("specify error.ratio argument"), smoothness = 0.5, t.type = "ar1", t.par = .5, D = NULL, T = NULL, reml = TRUE, lower = NULL, upper = NULL, control = list(trace = TRUE), optimizer="nlminb")
maxlik.cov.st(X, y, coords, time, sp.type = "exponential", range.par = stop("specify range.par argument"), error.ratio = stop("specify error.ratio argument"), smoothness = 0.5, t.type = "ar1", t.par = .5, D = NULL, T = NULL, reml = TRUE, lower = NULL, upper = NULL, control = list(trace = TRUE), optimizer="nlminb")
X |
A numeric matrix of size |
y |
A vector of length |
coords |
A numeric matrix of size |
time |
A numeric vector of length n containing the time at which the responses were observed. |
sp.type |
A character vector specifying the spatial covariance type. Valid types are currently exponential, gaussian, matern, and spherical. |
range.par |
An initial guess for the spatial dependence parameter. |
error.ratio |
A value non-negative value indicating the ratio |
smoothness |
A positive number indicating the variance of the error term. |
t.type |
A character vector indicating the spatial covariance type. Only |
t.par |
A value specifying the temporal dependence parameter of the ar1 process. |
D |
The Euclidean distance matrix for the coords matrix. Must be of size |
T |
The Euclidean distance matrix for the time matrix. Must be of size |
reml |
A boolean value indicating whether restricted maximum likelihood estimation should be used. Defaults to TRUE. |
lower |
A vector giving lower bounds for the covariance parameters |
upper |
A vector giving upper bounds for the covariance parameters |
control |
A list giving tuning parameters for the |
optimizer |
A vector describing the optimization function to use for the optimization. Currently, only |
When doing the numerical optimization, the covariance function is reparameterized slightly to speedup computation.
Specifically, the variance parameter for the process of interest,sp.par[1]
, is profiled out,
and the error.var
parameter is parameterized as sp.par[1] * error.ratio
, where error.ratio = error.var/sp.par[1]
.
Returns a list with the following elements:
sp.type |
The covariance form used. |
sp.par |
A vector containing the estimated variance of the hidden process and the spatial dependence. |
error.var |
The estimated error variance. |
smoothness |
The smoothness of the matern covariance function. |
par |
The final values of the optimization parameters. Note that these will not necessarily match |
convergence |
Convergence message from |
message |
Message from |
iterations |
Number of iterations for optimization to converge. |
evaluations |
Evaluations from |
Joshua French
cov.st
#Generate locations and observed times coords <- matrix(rnorm(40), ncol = 2) time <- rep(1:2, each = 10) #Calculate distance matrix for time vector T <- dist1(matrix(time)) #create design matrix X <- cbind(1, coords) #create mean for observed data to be generated mu <- X %*% c(1, 2, 3) #generate covariance matrix for spatio-temporal data V <- exp(-dist1(coords)) * .25^T #generate observe data y <- rmvnorm(mu = mu, V = V) maxlik.cov.st(X = X, y = y, coords = coords, time = time, sp.type = "exponential", range.par = 1, error.ratio = 0, t.type = "ar1", t.par = .5, reml = TRUE)
#Generate locations and observed times coords <- matrix(rnorm(40), ncol = 2) time <- rep(1:2, each = 10) #Calculate distance matrix for time vector T <- dist1(matrix(time)) #create design matrix X <- cbind(1, coords) #create mean for observed data to be generated mu <- X %*% c(1, 2, 3) #generate covariance matrix for spatio-temporal data V <- exp(-dist1(coords)) * .25^T #generate observe data y <- rmvnorm(mu = mu, V = V) maxlik.cov.st(X = X, y = y, coords = coords, time = time, sp.type = "exponential", range.par = 1, error.ratio = 0, t.type = "ar1", t.par = .5, reml = TRUE)
Plot contour lines from list produced by contourLines
function.
## S3 method for class 'contourLines' plot(x, begin=1, end = length(x), add = FALSE, ...)
## S3 method for class 'contourLines' plot(x, begin=1, end = length(x), add = FALSE, ...)
x |
The list of contour lines (created by |
begin |
Beginning position in list of contour lines you want to plot. |
end |
Ending position in list of contour lines you want to plot. |
add |
A boolean value indicating whether the contour lines should be added to an existing plot (add = TRUE) or should be plotted on a new plot (add = FALSE). |
... |
Additional arguments that will be passed to the |
This function does not return anything; it only creates a new plot or modifies an existing plot.
Joshua French
data(volcano) x <- 10*1:nrow(volcano) y <- 10*1:ncol(volcano) cL <- contourLines(x, y, volcano) plot.contourLines(cL)
data(volcano) x <- 10*1:nrow(volcano) y <- 10*1:ncol(volcano) cL <- contourLines(x, y, volcano) plot.contourLines(cL)
Generates realizations from a multivariate normal distribution conditional on observed data vector
rcondnorm(nsim = 1, y, mu, mup, V, Vp, Vop, method = "eigen")
rcondnorm(nsim = 1, y, mu, mup, V, Vp, Vop, method = "eigen")
nsim |
An integer indicating the number of realizations from the distribution. |
y |
A vector of length |
mu |
The mean vector of the observed data. Should be a vector of length |
mup |
The mean vector of the responses to be generated. Should be a vector of length |
V |
The covariance matrix of the observed data. The matrix should be symmetric and positive definite. The size must be |
Vp |
The covariance matrix of the responses to be generated. The matrix should be symmetric and positive definite. The size must be |
Vop |
The cross-covariance matrix between the observed data and the responses to be generated. The size must be |
method |
The method for performing a decomposition of the covariance matrix. Possible values are "eigen", "chol", and "svd", Eigen value decomposition, Cholesky decomposition, or Singular Value Decomposoition, respectively. |
An matrix containing the
nsim
realizations of the conditional normal distribution. Each column of the matrix represents a realization of the multivariate normal distribution.
Joshua French
rmvnorm
n <- 100 np <- 100 mu <- rep(1, 100) mup <- rep(2, 100) coords <- matrix(runif(2 * n), ncol = 2) pcoords <- matrix(runif(2 * np), ncol = 2) myV <- cov.sp(coords, sp.type = "exponential", c(1, 2), error.var = 1, pcoords = pcoords) y <- rmvnorm(1, mu = mu, V = myV$V) rcondnorm(3, y = y, mu = mu, mup = mup, V = myV$V, Vp = myV$Vp, Vop = myV$Vop, method = "chol")
n <- 100 np <- 100 mu <- rep(1, 100) mup <- rep(2, 100) coords <- matrix(runif(2 * n), ncol = 2) pcoords <- matrix(runif(2 * np), ncol = 2) myV <- cov.sp(coords, sp.type = "exponential", c(1, 2), error.var = 1, pcoords = pcoords) y <- rmvnorm(1, mu = mu, V = myV$V) rcondnorm(3, y = y, mu = mu, mup = mup, V = myV$V, Vp = myV$Vp, Vop = myV$Vop, method = "chol")
Generates realizations from a multivariate normal distribution.
rmvnorm(nsim = 1, mu, V, method = "eigen")
rmvnorm(nsim = 1, mu, V, method = "eigen")
nsim |
An integer indicating the number of realizations from the distribution. |
mu |
A vector of length n containing the mean values of the multivariate normal distribution. |
V |
The covariance matrix of the multivariate normal distribution. The matrix should be symmetric and positive definite. The size must be |
method |
The method for performing a decomposition of the covariance matrix. Possible values are "eigen", "chol", and "svd", Eigen value decomposition, Cholesky decomposition, or Singular Value Decomposoition, respectively. |
An matrix containing the
nsim
realizations of the multivariate normal distribution. Each column of the matrix represents a realization of the multivariate normal distribution.
Joshua French
rmvnorm
n <- 20 mu <- 1:n V <- exp(-dist1(matrix(rnorm(n)))) rmvnorm(nsim = 100, mu = mu, V = V, method = "eigen")
n <- 20 mu <- 1:n V <- exp(-dist1(matrix(rnorm(n)))) rmvnorm(nsim = 100, mu = mu, V = V, method = "eigen")
Calculates a spatial covariance using a (Euclidean) distance matrix D
. Not intended to be used directly by user (though it may be helpful to some). It is strongly recommended that you use the cov.sp
function. No argument or error checking is provided for this function.
simple.cov.sp(D, sp.type, sp.par, error.var, smoothness, finescale.var)
simple.cov.sp(D, sp.type, sp.par, error.var, smoothness, finescale.var)
D |
A distance matrix between locations |
sp.type |
A character vector specifying the spatial covariance type. Valid types are currently exponential, gaussian, matern, and spherical. |
sp.par |
A vector of length 2 specifying the scale and dependence of the covariance function. The first element refers to the variance of the hidden process (sometimes this is called the partial sill) while the second elements determines the strength of dependence between locations. |
error.var |
A non-negative number indicating the variance of the error term. |
smoothness |
A positive number indicating the variance of the error term. |
finescale.var |
A non-negative positive number indicating the finescale variability. The is also called the microscale variance |
Returns a covariance matrix.
Joshua French
~ cov.sp
coords <- matrix(rnorm(30), ncol = 3) D <- dist1(coords) simple.cov.sp(D = D, sp.type = "exponential", sp.par = c(2, 1), error.var = 1, smoothness = 0.5, finescale.var = 0)
coords <- matrix(rnorm(30), ncol = 3) D <- dist1(coords) simple.cov.sp(D = D, sp.type = "exponential", sp.par = c(2, 1), error.var = 1, smoothness = 0.5, finescale.var = 0)
Calculates a temporal covariance using a (Euclidean) distance matrix T
. Not intended to be used directly by user (though it may be helpful to some). It is used in the covets
function. No argument or error checking is provided for this function.
simple.cov.time(T, t.type, t.par)
simple.cov.time(T, t.type, t.par)
T |
A distance matrix. |
t.type |
A character vector specifying the temporal covariance type. Only "ar1" is currently implemented. |
t.par |
A vector of length 1 specifying the strength of dependence of the covariance function. |
Returns a covariance matrix.
Joshua French
~ cov.st
T <- dist1(matrix(1:10)) simple.cov.time(T = T, t.type = "ar1", t.par = .5)
T <- dist1(matrix(1:10)) simple.cov.time(T = T, t.type = "ar1", t.par = .5)
The function spLMPredictJoint
collects posterior predictive samples
for a set of new locations given a spLM
object from the spBayes
package.
spLMPredictJoint(sp.obj, pred.coords, pred.covars, start = 1, end = nrow(sp.obj$p.theta.samples), thin = 1, verbose = TRUE, n.report = 100, noisy = FALSE, method = "eigen")
spLMPredictJoint(sp.obj, pred.coords, pred.covars, start = 1, end = nrow(sp.obj$p.theta.samples), thin = 1, verbose = TRUE, n.report = 100, noisy = FALSE, method = "eigen")
sp.obj |
An |
pred.coords |
An |
pred.covars |
An |
start |
Specifies the first sample included in the composition sampling. |
end |
Specifies the last sample included in the composition.
The default is to use all posterior samples in |
thin |
A sample thinning factor. The default of 1 considers all
samples between |
verbose |
If |
n.report |
The interval to report sampling progress. |
noisy |
If |
method |
Method used to decompose covariance matrix. Options are "chol", "eigen", and "svd" for the Cholesky, Eigen, and singular value decomposition approaches, respectively. |
This function samples from the joint posterior predictive distribution of a Bayesian spatial linear model. Specifically, it is intended to be similar to the spPredict
function in the spBayes
except that it samples from the joint distribution instead of the marginal distribution. However, it will only work for spLM
objects and should have the same limitations as the spLM
and spPredict
functions. Note that the spRecover
function is called internally to recover the posterior samples form the posterior distribution of the spatial model.
The function returns a matrix of posterior predictive samples, where
B
is the number of posterior samples. The class is jointPredicitveSample
.
Joshua French
spLM, spPredict, spRecover
# Set parameters n <- 100 np <- 12 n.samples <- 10 burnin.start <- .5 * n.samples + 1 sigmasq <- 1 tausq <- 0.0 phi <- 1 cov.model <- "exponential" n.report <- 5 # Generate coordinates coords <- matrix(runif(2 * n), ncol = 2); pcoords <- as.matrix(expand.grid(seq(0, 1, len = 12), seq(0, 1, len = 12))) # Construct design matrices X <- as.matrix(cbind(1, coords)) Xp <- cbind(1, pcoords) # Specify priors starting <- list("phi" = phi, "sigma.sq"= sigmasq, "tau.sq" = tausq) tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1) priors.1 <- list("beta.Norm"=list(c(1, 2, 1), diag(100, 3)), "phi.Unif"=c(0.00001, 10), "sigma.sq.IG"=c(1, 1)) # Generate data B <- rnorm(3, c(1, 2, 1), sd = 10) phi <- runif(1, 0, 10) sigmasq <- 1/rgamma(1, 1, 1) V <- simple.cov.sp(D = dist1(coords), cov.model, c(sigmasq, 1/phi), error.var = tausq, smoothness = nu, finescale.var = 0) y <- X %*% B + rmvnorm(1, rep(0, n), V) + rnorm(n, 0, sqrt(tausq)) # Create spLM object library(spBayes) m1 <- spBayes::spLM(y ~ X - 1, coords = coords, starting = starting, tuning = tuning, priors = priors.1, cov.model = cov.model, n.samples = n.samples, verbose = FALSE, n.report = n.report) # Sample from joint posterior predictive distribution y1 <- spLMPredictJoint(m1, pred.coords = pcoords, pred.covars = Xp, start = burnin.start, verbose = FALSE, method = "chol")
# Set parameters n <- 100 np <- 12 n.samples <- 10 burnin.start <- .5 * n.samples + 1 sigmasq <- 1 tausq <- 0.0 phi <- 1 cov.model <- "exponential" n.report <- 5 # Generate coordinates coords <- matrix(runif(2 * n), ncol = 2); pcoords <- as.matrix(expand.grid(seq(0, 1, len = 12), seq(0, 1, len = 12))) # Construct design matrices X <- as.matrix(cbind(1, coords)) Xp <- cbind(1, pcoords) # Specify priors starting <- list("phi" = phi, "sigma.sq"= sigmasq, "tau.sq" = tausq) tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1) priors.1 <- list("beta.Norm"=list(c(1, 2, 1), diag(100, 3)), "phi.Unif"=c(0.00001, 10), "sigma.sq.IG"=c(1, 1)) # Generate data B <- rnorm(3, c(1, 2, 1), sd = 10) phi <- runif(1, 0, 10) sigmasq <- 1/rgamma(1, 1, 1) V <- simple.cov.sp(D = dist1(coords), cov.model, c(sigmasq, 1/phi), error.var = tausq, smoothness = nu, finescale.var = 0) y <- X %*% B + rmvnorm(1, rep(0, n), V) + rnorm(n, 0, sqrt(tausq)) # Create spLM object library(spBayes) m1 <- spBayes::spLM(y ~ X - 1, coords = coords, starting = starting, tuning = tuning, priors = priors.1, cov.model = cov.model, n.samples = n.samples, verbose = FALSE, n.report = n.report) # Sample from joint posterior predictive distribution y1 <- spLMPredictJoint(m1, pred.coords = pcoords, pred.covars = Xp, start = burnin.start, verbose = FALSE, method = "chol")
A list containing X
, a 50 x 3 design matrix, y
, a vector of length 50 of observed responses, V
, a 50 x 50 covariance matrix for the observed data, Xp
, a 121 x 3 design matrix for the predicted responses, Vp
, the 121 x 121 covariance matrix of the predicted responses, Vop
, the 50 x 121 covariance matrix between the observed responses and the predicted responses, coords
, a 50 x 2 matrix containing the sites of the 50 observed responses, and pcoords
, a 121 x 2 matrix containing the 121 sites for the predicted responses (a 11 x 11 regular grid over the domain [0, 1]x[0, 1]).
data(toydata)
data(toydata)
Joshua French
data(toydata)
data(toydata)