predict.NNGP | R Documentation |
NNGP
models.The function predict
collects posterior predictive samples
for a set of new locations given an object of
class NNGP
.
## S3 method for class 'NNGP' predict(object, X.0, coords.0, sub.sample, n.omp.threads = 1, verbose=TRUE, n.report=100, ...)
object |
an object of class |
X.0 |
the design matrix for prediction locations. An
intercept should be provided in the first column if one is specified
in |
coords.0 |
the spatial coordinates corresponding to
|
sub.sample |
an optional list that specifies the samples to included in
the composition sampling a non-Conjugate model. Valid tags are |
n.omp.threads |
a positive integer indicating
the number of threads to use for SMP parallel processing. The package must
be compiled for OpenMP support. For most Intel-based machines, we recommend setting
|
verbose |
if |
n.report |
the interval to report sampling progress. |
... |
currently no additional arguments. |
An object of class predict.NNGP
which is a list comprising:
p.y.0 |
a matrix that holds the response variable posterior
predictive samples where rows are locations corresponding to
|
p.w.0 |
a matrix that holds the random effect posterior
predictive samples where rows are locations corresponding to
|
run.time |
execution time reported using |
Andrew O. Finley finleya@msu.edu,
Abhirup Datta abhidatta@jhu.edu,
Sudipto Banerjee sudipto@ucla.edu
Datta, A., S. Banerjee, A.O. Finley, and A.E. Gelfand. (2016) Hierarchical Nearest-Neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, doi: 10.1080/01621459.2015.1044091.
Finley, A.O., A. Datta, B.D. Cook, D.C. Morton, H.E. Andersen, and S. Banerjee. (2019) Efficient algorithms for Bayesian Nearest Neighbor Gaussian Processes. Jurnal of Computational and Graphical Statistics, doi: 10.1080/10618600.2018.1537924.
rmvn <- function(n, mu=0, V = matrix(1)){ p <- length(mu) if(any(is.na(match(dim(V),p)))) stop("Dimension problem!") D <- chol(V) t(matrix(rnorm(n*p), ncol=p)%*%D + rep(mu,rep(n,p))) } ##Make some data set.seed(1) n <- 100 coords <- cbind(runif(n,0,1), runif(n,0,1)) x <- cbind(1, rnorm(n)) B <- as.matrix(c(1,5)) sigma.sq <- 5 tau.sq <- 1 phi <- 3/0.5 D <- as.matrix(dist(coords)) R <- exp(-phi*D) w <- rmvn(1, rep(0,n), sigma.sq*R) y <- rnorm(n, x%*%B + w, sqrt(tau.sq)) ho <- sample(1:n, 50) y.ho <- y[ho] x.ho <- x[ho,,drop=FALSE] w.ho <- w[ho] coords.ho <- coords[ho,] y <- y[-ho] x <- x[-ho,,drop=FALSE] w <- w[-ho,,drop=FALSE] coords <- coords[-ho,] ##Fit a Response, Latent, and Conjugate NNGP model n.samples <- 500 starting <- list("phi"=phi, "sigma.sq"=5, "tau.sq"=1) tuning <- list("phi"=0.5, "sigma.sq"=0.5, "tau.sq"=0.5) priors <- list("phi.Unif"=c(3/1, 3/0.01), "sigma.sq.IG"=c(2, 5), "tau.sq.IG"=c(2, 1)) cov.model <- "exponential" n.report <- 500 ##Latent m.s <- spNNGP(y~x-1, coords=coords, starting=starting, method="latent", n.neighbors=10, tuning=tuning, priors=priors, cov.model=cov.model, n.samples=n.samples, n.omp.threads=1, n.report=n.report) p.s <- predict(m.s, X.0 = x.ho, coords.0 = coords.ho, n.omp.threads=1) plot(apply(p.s$p.w.0, 1, mean), w.ho) plot(apply(p.s$p.y.0, 1, mean), y.ho) ##Response m.r <- spNNGP(y~x-1, coords=coords, starting=starting, method="response", n.neighbors=10, tuning=tuning, priors=priors, cov.model=cov.model, n.samples=n.samples, n.omp.threads=1, n.report=n.report) p.r <- predict(m.r, X.0 = x.ho, coords.0 = coords.ho, n.omp.threads=1) points(apply(p.r$p.y.0, 1, mean), y.ho, pch=19, col="blue") ##Conjugate theta.alpha <- c(phi, tau.sq/sigma.sq) names(theta.alpha) <- c("phi", "alpha") m.c <- spConjNNGP(y~x-1, coords=coords, n.neighbors=10, theta.alpha=theta.alpha, sigma.sq.IG=c(2, sigma.sq), cov.model=cov.model, n.omp.threads=1) p.c <- predict(m.c, X.0 = x.ho, coords.0 = coords.ho, n.omp.threads=1) points(p.c$y.0.hat, y.ho, pch=19, col="orange")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.