geoBayes (version 0.3.0)

bf2new: Compute the Bayes factors at new points

Description

Compute the Bayes factors.

Usage

bf2new(bf1obj, linkp, phi, omg, kappa, useCV = TRUE)

Arguments

bf1obj

Output from the function bf1skel which contains the Bayes factors and importance sampling weights.

linkp,phi,omg,kappa

Optional scalar or vector or NULL. If scalar or vector, the Bayes factors are calculated at those values with respect to the reference model used in bf1skel. If missing or NULL the

useCV

Whether to use control variates for finer corrections.

Value

• An array of size length(linkp) * length(phi) * length(omg) * length(kappa) containing the Bayes factors for each combination of the parameters.

Details

Computes the Bayes factors using the importance weights at the new points. The new points are taken from the grid derived by expanding the parameter values inputted. The arguments linkp phi omg kappa correspond to the link function, spatial range, relative nugget, and correlation function parameters respectively.

References

Doss, H. (2010). Estimation of large families of Bayes factors from Markov chain output. Statistica Sinica, 20(2), 537.

Roy, V., Evangelou, E., and Zhu, Z. (2014). Efficient estimation and prediction for the Bayesian spatial generalized linear mixed model with flexible link functions. Technical report, Iowa State University.

Examples

data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
linkp <- "probit"
### Skeleton points
philist <- c(100, 140, 180)
omglist <- c(.5, 1)
parlist <- expand.grid(phi=philist, linkp=linkp, omg=omglist, kappa = kappa)
### MCMC sizes
Nout <- 100
Nthin <- 1
Nbi <- 0
### Take MCMC samples
runs <- list()
for (i in 1:NROW(parlist)) {
  runs[[i]] <- mcsglmm(Infected ~ 1, 'binomial', rhizoctonia, weights = Total,
                       atsample = ~ Xcoord + Ycoord,
                       Nout = Nout, Nthin = Nthin, Nbi = Nbi,
                       betm0 = betm0, betQ0 = betQ0,
                       ssqdf = ssqdf, ssqsc = ssqsc,
                       phistart = parlist$phi[i], omgstart = parlist$omg[i],
                       linkp = parlist$linkp[i], kappa = parlist$kappa[i],
                       corrfcn = corrf, phisc = 0, omgsc = 0)
}
bf <- bf1skel(runs)
bfall <- bf2new(bf, phi = seq(100, 200, 10), omg = seq(0, 2, .2))
plotbf2(bfall, c("phi", "omg"))