spDiag | R Documentation |
The function spDiag
calculates DIC, GP, GRS, and associated
statistics given a spLM
, spMvLM
,
spGLM
, spMvGLM
, spMvGLM
, or
spSVC
object.
spDiag(sp.obj, start=1, end, thin=1, verbose=TRUE, n.report=100, ...)
sp.obj |
an object returned by |
start |
specifies the first sample included in the computation. The |
end |
specifies the last sample included in the computation.
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 progress. |
... |
currently no additional arguments. |
A list with some of the following tags:
DIC |
a matrix holding DIC and associated statistics, see Banerjee et al. (2004) for details. |
GP |
a matrix holding GP and associated statistics, see Gelfand and Ghosh (1998) for details. |
GRS |
a scoring rule, see Equation 27 in Gneiting and Raftery (2007) for details. |
Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee sudipto@ucla.edu
Banerjee, S., Carlin, B.P., and Gelfand, A.E. (2004). Hierarchical modeling and analysis for spatial data. Chapman and Hall/CRC Press, Boca Raton,Fla.
Finley, A.O. and S. Banerjee (2019) Efficient implementation of spatially-varying coefficients models.
Gelfand A.E. and Ghosh, S.K. (1998). Model choice: a minimum posterior predictive loss approach. Biometrika. 85:1-11.
Gneiting, T. and Raftery, A.E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association. 102:359-378.
## Not run: 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))) } set.seed(1) n <- 100 coords <- cbind(runif(n,0,1), runif(n,0,1)) X <- as.matrix(cbind(1, rnorm(n))) B <- as.matrix(c(1,5)) p <- length(B) sigma.sq <- 2 tau.sq <- 0.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)) n.samples <- 1000 starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1) tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1) ##too restrictive of prior on beta priors.1 <- list("beta.Norm"=list(rep(0,p), diag(1,p)), "phi.Unif"=c(3/1, 3/0.1), "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1)) ##more reasonable prior for beta priors.2 <- list("beta.Norm"=list(rep(0,p), diag(1000,p)), "phi.Unif"=c(3/1, 3/0.1), "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1)) cov.model <- "exponential" n.report <- 500 verbose <- TRUE m.1 <- spLM(y~X-1, coords=coords, starting=starting, tuning=tuning, priors=priors.1, cov.model=cov.model, n.samples=n.samples, verbose=verbose, n.report=n.report) m.2 <- spLM(y~X-1, coords=coords, starting=starting, tuning=tuning, priors=priors.2, cov.model=cov.model, n.samples=n.samples, verbose=verbose, n.report=n.report) ##non-spatial model m.3 <- spLM(y~X-1, n.samples=n.samples, verbose=verbose, n.report=n.report) burn.in <- 0.5*n.samples ##recover beta and spatial random effects m.1 <- spRecover(m.1, start=burn.in, verbose=FALSE) m.2 <- spRecover(m.2, start=burn.in, verbose=FALSE) ##lower is better for DIC, GPD, and GRS print(spDiag(m.1)) print(spDiag(m.2)) print(spDiag(m.3)) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.