Compute the Bayes factors at new points

Share:

Description

Compute the Bayes factors.

Usage

1
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 then the unique values from the MCMC chains that were inputted in bf1skel will be used.

useCV

Whether to use control variates for finer corrections.

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.

Value

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

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. (2015). Efficient estimation and prediction for the Bayesian spatial generalized linear mixed model with flexible link functions. Biometrics. http://dx.doi.org/10.1111/biom.12371

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## Not run: 
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"))

## End(Not run)