Model fit diagnostics

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 spLM, spMvLM, spGLM, spMvGLM. For spSVC, sp.obj is an object from spRecover.

start

specifies the first sample included in the computation. The start, end, and thin arguments only apply to spGLM or spMvGLM objects. Sub-sampling for spLM and spMvLM is controlled using spRecover which must be called prior to spDiag.

end

specifies the last sample included in the computation. The default is to use all posterior samples in sp.obj. See start argument description.

thin

a sample thinning factor. The default of 1 considers all samples between start and end. For example, if thin = 10 then 1 in 10 samples are considered between start and end.

verbose

if TRUE calculation progress is printed to the screen; otherwise, nothing is printed to the screen.

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.

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.

Author

Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee sudipto@ucla.edu

Examples

if (FALSE) { # \dontrun{
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))
} # }