spDiag {spBayes} | R Documentation |
Model fit diagnostics
Description
The function spDiag
calculates DIC, GP, GRS, and associated
statistics given a spLM
, spMvLM
,
spGLM
, spMvGLM
, spMvGLM
, or
spSVC
object.
Usage
spDiag(sp.obj, start=1, end, thin=1, verbose=TRUE, n.report=100, ...)
Arguments
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. |
Value
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. |
Author(s)
Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee sudipto@ucla.edu
References
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.
Examples
## 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)