WARNING: RUNNING THIS IS VERY TIME-CONSUMING

> ##
> ## 1. Simulating data
> ##
> sim <- grf(grid = expand.grid(x = seq(1.1, 7.1, l = 7), y = seq(1.1, 7.1, l = 7)), cov.pars = c(0.1, 2))
grf: process with 1 covariance structure(s)
grf: nugget effect is: tausq= 0
grf: covariance model 1 is: exponential(sigmasq=0.1, phi=2)
grf: decomposition algorithm used is: cholesky
grf: simulation(s) on a grid provided by the user
grf: End of simulation procedure. Number of realizations: 1
> sim$units.m <- rep(4,49)
> sim$data <- rbinom(49, size = rep(4,49), prob = exp(sim$data)/(1+exp(sim$data)))
>
> ## 2. Visualising the data
> ##
> plot(sim$coords, type = "n")
> text(sim$coords[,1], sim$coords[,2], format(sim$data), cex=1.5)

> ## 3. Setting input options
>
> sim.pr <- prior.glm.control(phi.discrete = seq(0.05, 3, l=60))
> sim.mcmc <- mcmc.control(S.scale = 0.1, phi.scale = 0.04, burn.in=10000)
> grid <- expand.grid(x = seq(1, 7, l = 51), y = seq(1, 7, l = 51))
>
> run.sim <- binom.krige.bayes(sim, locations = grid, prior = sim.pr, mcmc.input = sim.mcmc)
binom.krige.bayes: analysis assuming constant mean
burn-in = 10000 is finished; Acc.-rate = 0.67 ; Acc-rate-phi = 0.72
iter. numb. 11000 ; Acc.-rate = 0.77 ; Acc-rate-phi = 0.69
iter. numb. 12000 ; Acc.-rate = 0.69 ; Acc-rate-phi = 0.67
iter. numb. 13000 ; Acc.-rate = 0.71 ; Acc-rate-phi = 0.65
iter. numb. 14000 ; Acc.-rate = 0.24 ; Acc-rate-phi = 0.59
iter. numb. 15000 ; Acc.-rate = 0.67 ; Acc-rate-phi = 0.70
iter. numb. 16000 ; Acc.-rate = 0.78 ; Acc-rate-phi = 0.76
iter. numb. 17000 ; Acc.-rate = 0.75 ; Acc-rate-phi = 0.73
iter. numb. 18000 ; Acc.-rate = 0.71 ; Acc-rate-phi = 0.67
iter. numb. 19000 ; Acc.-rate = 0.71 ; Acc-rate-phi = 0.67
iter. numb. 20000 ; Acc.-rate = 0.67 ; Acc-rate-phi = 0.64
MCMC performed: n.iter. = 10000 ; thinning = 10 ; burn.in = 10000
binom.krige.bayes: Prediction performed
binom.krige.bayes: done!
>
> ## inspecting output
> names(run.sim)
[1] "posterior" "predictive" "model" "prior" "mcmc.input"
[6] ".Random.seed" "call"
> names(run.sim$posterior)
[1] "phi" "beta" "sigmasq" "simulations"
> names(run.sim$posterior$phi)
[1] "sample" "mean" "var"
>
> par(mfrow=c(2,2), mar=c(2.3,2.5,.5,.7), mgp=c(1.5,.6,0), cex=0.6)
> plot(run.sim$posterior$phi$sample,type="l")
> acf(run.sim$posterior$phi$sample)
> plot(run.sim$posterior$sim[1,],type="l")
> acf(run.sim$posterior$sim[1,])

> ## Predictions
> ##
>
> names(run.sim$predictive)
[1] "simulations" "median" "uncertainty" "quantiles"
>
> ##
> ## Plotting predictions
> ##
> image.kriging(locations = grid, values = run.sim$pred$median )