The geoBayes package

Description

Analysis of geostatistical data using Bayes and Empirical Bayes methods.

Details

This package provides functions to fit geostatistical data. The data can be continuous, binary or count data and the models implemented are flexible. Conjugate priors are assumed on some parameters while inference on the other parameters can be done through a full Bayesian analysis of by empirical Bayes methods.

Some demonstration examples are provided. Type demo(package = "geoBayes") to examine them.

Author(s)

Evangelos Evangelou <e.evangelou@maths.bath.ac.uk> and Vivekananda Roy <vroy@iastate.edu>

References

Roy, V., Evangelou, E. and Zhu, Z. (2014). Empirical Bayes methods for the transformed Gaussian random fields model with additive measurement errors. In Upadhyay, S. K., Singh, U., Dey, D. K., and Loganathan, A., editors, Current Trends in Bayesian Methodology with Applications, Boca Raton, FL, USA, CRC Press.

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, 72(1), 289-298.

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Roy, V., & Evangelou, E. (2024). Selection of proposal distributions for multiple importance sampling. Statistica Sinica, 34, 27-46.

Examples

## Not run: 
demo(package = "geoBayes")
demo(rhizoctonia1, package = "geoBayes")
demo(rhizoctonia1, package = "geoBayes")

## End(Not run)

Approximate log-likelihood calculation

Description

Calculate the likelihood approximation at different parameter values. This function is useful for choosing the skeleton set.

Plot likelihood approximation.

Usage

alik_cutoff(likopt, par_vals, likthreshold)

alik_plot(alikobj)

Arguments

Details

The input par_vals is meant to contain vector of parameter values for each parameter. For each element in par_vals, the other parameters are set equal to the maximisers given in likopt and the approximate likelihood is computed. The cuttoff is calculated using linear interpolation provided by approx.

The plot can be used to visualise the Laplace approximation to the likelihood provided by the function alik_cutoff.

Value

A list with the log-likelihood approximation and cutoff values.

Draws a plot.

References

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Log-likelihood approximation

Description

Log-likelihood approximation.

Usage

alik_inla(
  par_vals,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  np,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  tsqdf,
  tsqsc,
  dispersion = 1,
  longlat = FALSE
)

Arguments

Details

Computes and approximation to the log-likelihood for the given parameters using integrated nested Laplace approximations.

Value

A list with components

• par_vals A data frame of the parameter values.

• aloglik The approximate log-likelihood at thos parameter values.

References

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Log-likelihood maximisation

Description

Approximate log-likelihood maximisation

Usage

alik_optim(
  paroptim,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  np,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  dispersion = 1,
  longlat = FALSE,
  control = list()
)

Arguments

Details

Uses the "L-BFGS-B" method of the function optim to maximise the log-likelihood for the parameters linkp, phi, omg, kappa.

Value

The output from the function optim. The "value" element is the log-likelihood, not the negative log-likelihood.

References

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Computation of Bayes factors at the skeleton points

Description

Function to compute the Bayes factors from MCMC samples.

Usage

bf1skel(
  runs,
  bfsize1 = 0.8,
  method = c("RL", "MW"),
  reference = 1,
  transf = c("no", "mu", "wo")
)

Arguments

Details

Computes the Bayes factors using method with respect to reference.

Value

A list with components

• logbf A vector containing logarithm of the Bayes factors.

• logLik1 logLik2 Matrices with the values of the log-likelihood computed from the samples for each model at the first and second stages.

• isweights A vector with the importance sampling weights for computing the Bayes factors at new points that will be used at the second stage. Used internally in bf2new and bf2optim.

• controlvar A matrix with the control variates computed at the samples that will be used in the second stage.

• sample2 The MCMC sample for mu or z that will be used in the second stage. Used internally in bf2new and bf2optim.

• N1, N2 Vectors containing the sample sizes used in the first and second stages.

• distmat Matrix of distances between locations.

• betm0, betQ0, ssqdf, ssqsc, tsqdf, tsqsc, dispersion, response, weights, modelmatrix, locations, family, corrfcn, transf Model parameters used internally in. bf2new and bf2optim.

• pnts A list containing the skeleton points. Used internally in bf2new and bf2optim.

References

Geyer, C. J. (1994). Estimating normalizing constants and reweighting mixtures. Technical report, University of Minnesota.

Meng, X. L., & Wong, W. H. (1996). Simulating ratios of normalizing constants via a simple identity: A theoretical exploration. Statistica Sinica, 6, 831-860.

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, 72(1), 289-298.

Examples

## Not run: 
data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
family <- "binomial.probit"
### Skeleton points
philist <- c(100, 140, 180)
omglist <- c(.5, 1)
parlist <- expand.grid(linkp=0, phi=philist, 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, family, rhizoctonia, weights = Total,
                       atsample = ~ Xcoord + Ycoord,
                       Nout = Nout, Nthin = Nthin, Nbi = Nbi,
                       betm0 = betm0, betQ0 = betQ0,
                       ssqdf = ssqdf, ssqsc = ssqsc,
                       phi = parlist$phi[i], omg = parlist$omg[i],
                       linkp = parlist$linkp[i], kappa = parlist$kappa[i],
                       corrfcn = corrf,
                       corrtuning=list(phi = 0, omg = 0, kappa = 0))
}
bf <- bf1skel(runs)
bf$logbf

## End(Not run)

Compute the Bayes factors at new points

Description

Compute the Bayes factors.

Usage

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

Arguments

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, 72(1), 289-298.

Examples

## Not run: 
data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
family <- "binomial.probit"
### Skeleton points
philist <- c(100, 140, 180)
omglist <- c(.5, 1)
parlist <- expand.grid(linkp=0, phi=philist, 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, family, rhizoctonia, weights = Total,
                       atsample = ~ Xcoord + Ycoord,
                       Nout = Nout, Nthin = Nthin, Nbi = Nbi,
                       betm0 = betm0, betQ0 = betQ0,
                       ssqdf = ssqdf, ssqsc = ssqsc,
                       phi = parlist$phi[i], omg = parlist$omg[i],
                       linkp = parlist$linkp[i], kappa = parlist$kappa[i],
                       corrfcn = corrf,
                       corrtuning=list(phi = 0, omg = 0, kappa = 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)

Empirical Bayes estimator

Description

Estimation by empirical Bayes.

Usage

bf2optim(bf1obj, paroptim, useCV = TRUE, control = list())

Arguments

Details

This function is a wrap around bf2new using the "L-BFGS-B" method of the function optim to estimate the parameters.

Value

The output from the function optim. The "value" element is the log-Bayes factor, not the negative log-Bayes factor.

Examples

## Not run: 
data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
family <- "binomial.probit"
### Skeleton points
philist <- c(100, 140, 180)
omglist <- c(.5, 1)
parlist <- expand.grid(linkp=0, phi=philist, 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, family, rhizoctonia, weights = Total,
                       atsample = ~ Xcoord + Ycoord,
                       Nout = Nout*c(.8,.2), Nthin = Nthin, Nbi = Nbi,
                       betm0 = betm0, betQ0 = betQ0,
                       ssqdf = ssqdf, ssqsc = ssqsc,
                       phi = parlist$phi[i], omg = parlist$omg[i],
                       linkp = parlist$linkp[i], kappa = parlist$kappa[i],
                       corrfcn = corrf,
                       corrtuning=list(phi = 0, omg = 0, kappa = 0))
}
bf <- bf1skel(runs)
est <- bf2optim(bf, list(phi = c(100, 200), omg = c(0, 2)))
est

## End(Not run)

Empirical Bayes standard errors

Description

Standard errors for the empirical Bayes estimates of the parameters.

Usage

bf2se(mcrun, transf = c("no", "mu", "wo"))

Arguments

Value

The precision (inverse covariance) matrix.

References

Casella, G. (2001). Empirical Bayes Gibbs sampling. Biostatistics, 2(4), 485-500.

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Batch means, Bayes factors standard errors

Description

Compute the standard errors for the Bayes factors estimates.

Usage

bmbfse(
  pargrid,
  runs,
  bfsize1 = 0.8,
  nbatch1 = 0.5,
  nbatch2 = 0.5,
  S1method = c("RL", "MW"),
  bvmethod = c("Standard", "TukeyHanning", "Bartlett"),
  reference = 1,
  transf = c("no", "mu", "wo")
)

Arguments

Details

Uses the batch means method to compute the standard errors for Bayes factors.

Value

A list with components

• pargrid The inputted pargrid augmented with the computed log Bayes factors and standard errors.

• bfEstimate The estimates of the Bayes factors

• bfSigma The covariance matrix for the Bayes factors estimates.

References

Roy, V., Tan, A. and Flegal, J. (2018). Estimating standard errors for importance sampling estimators with multiple Markov chains, Statistica Sinica, 28 1079-1101.

Roy, V., & Evangelou, E. (2018). Selection of proposal distributions for generalized importance sampling estimators. arXiv preprint arXiv:1805.00829.

Spatial correlation used in the geoBayes package

Description

This hidden variable contains a choice of correlation functions that can be fit with this package. The function can be chosen in the corrfcn input of the relevant function. This variable cannot be overwritten.

Usage

.geoBayes_corrfcn

Format

An object of class data.frame with 5 rows and 4 columns.

Models used in the geoBayes package

Description

This hidden variable contains a choice of models that can be fit with this package. The model can be chosen in the family input of the relevant function. This variable cannot be overwritten.

Usage

.geoBayes_models

Format

An object of class data.frame with 12 rows and 7 columns.

Calculate the link function for exponential families

Description

Link function for the exponential family.

Usage

linkfcn(mu, linkp, family = "gaussian")

linkinv(z, linkp, family = "gaussian")

Arguments

Details

linkfcn maps the mean of the response variable mu to the linear predictor z. linkinv is its inverse.

For the Gaussian family, if the link parameter is positive, then the extended link is used, defined by

z = \frac{sign(\mu)|\mu|^\nu - 1}{\nu}

In the other case, the usual Box-Cox link is used.

For the Poisson and gamma families, if the link parameter is positive, then the link is defined by

z = \frac{sign(w) (e^{\nu |w|}-1)}{\nu}

where w = \log(\mu). In the other case, the usual Box-Cox link is used.

For the GEV binomial family, the link function is defined by

\mu = 1 - \exp\{-\max(0, 1 + \nu z)^{\frac{1}{\nu}}\}

for any real \nu. At \nu = 0 it reduces to the complementary log-log link.

The Wallace binomial family is a fast approximation to the robit family. It is defined as

\mu = \Phi(\mbox{sign}(z) c(\nu) \sqrt{\nu \log(1 + z^2/\nu)})

where c(\nu) = (8\nu+1)/(8\nu+3)

Value

A numeric array of the same dimension as the function's first argument.

References

Evangelou, E., & Roy, V. (2019). Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling. Spatial Statistics, 29, 289-315.

Examples

## Not run: 
mu <- seq(0.1, 0.9, 0.1)
linkfcn(mu, 7, "binomial")       # robit(7) link function
linkfcn(mu, , "binomial.logit")  # logit link function

mu <- seq(-3, 3, 1)
linkfcn(mu, 0.5, "gaussian")     # sqrt transformation
linkinv(linkfcn(mu, 0.5, "gaussian"), 0.5, "gaussian")
curve(linkfcn(x, 0.5, "gaussian"), -3, 3)

## End(Not run)

Convert to an mcmc object

Description

Convert to an mcmc object.

Usage

mcmcmake(...)

Arguments

Details

This function takes as input the one or more output(s) from function mcsglmm or mcstrga and returns an mcmc object or an mcmc.list object for coda. The function requires the coda package to be installed.

The spatial random field components are assigned the names z_* where * is a number beginning at 1. Similarly, the regressor coefficients are assigned the names beta_* if not unique, or simply beta if there is only one regressor. The names ssq, tsq, phi, omg correspond to the partial sill, measurement error variance, spatial range, and relative nugget parameters respectively.

Value

An mcmc object.

See Also

Functions such as plot.mcmc and summary.mcmc in the coda package. The function do.call can be used to pass arguments stored in a list.

Examples

## Not run: 
 ### Load the data
data(rhizoctonia)
rhiz <- na.omit(rhizoctonia)
rhiz$IR <- rhiz$Infected/rhiz$Total # Incidence rate of the
                              # rhizoctonia disease
 ### Define the model
corrf <- "spherical"
ssqdf <- 1
ssqsc <- 1
tsqdf <- 1
tsqsc <- 1
betm0 <- 0
betQ0 <- diag(.01, 2, 2)
phiprior <- c(200, 1, 1000, 100) # U(100, 300)
phisc <- 1
omgprior <- c(3, 1, 1000, 0) # U(0, 3)
omgsc <- 1.3
linkp <- 1
## MCMC parameters
Nout <- 100
Nbi <- 0
Nthin <- 1
 ### Run MCMC
sample <- mcstrga(Yield ~ IR, data = rhiz,
                  atsample = ~ Xcoord + Ycoord, corrf = corrf,
                  Nout = Nout, Nthin = Nthin,
                  Nbi = Nbi, betm0 = betm0, betQ0 = betQ0,
                  ssqdf = ssqdf, ssqsc = ssqsc,
                  tsqdf = tsqdf, tsqsc = tsqsc,
                  linkp = linkp,
                  corrprior = list(phi = phiprior, omg = omgprior), 
                  corrtuning = list(phi = phisc, omg = omgsc, kappa = 0),
                  test = FALSE)
mcsample <- mcmcmake(sample)
plot(mcsample[, c("phi", "omg", "beta_1", "beta_2", "ssq", "tsq")],
     density = FALSE)
summary(mcsample[, c("phi", "omg", "beta_1", "beta_2", "ssq", "tsq")])

## End(Not run)

MCMC samples from the Spatial GLMM

Description

Draw MCMC samples from the Spatial GLMM with known link function

Usage

mcsglmm(
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  linkp,
  phi,
  omg,
  kappa,
  Nout,
  Nthin = 1,
  Nbi = 0,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  corrpriors,
  corrtuning,
  dispersion = 1,
  longlat = FALSE,
  test = FALSE
)

Arguments

Details

The four-parameter prior for phi is defined by

\propto (\phi - \theta_4)^{\theta_2 -1} \exp\{-(\frac{\phi - \theta_4}{\theta_1})^{\theta_3}\}

for \phi > \theta_4. The prior for omg is similar. The prior parameters correspond to scale, shape, exponent, and location. See arXiv:1005.3274 for details of this distribution.

The GEV (Generalised Extreme Value) link is defined by

\mu = 1 - \exp\{-\max(0, 1 + \nu x)^{\frac{1}{\nu}}\}

for any real \nu. At \nu = 0 it reduces to the complementary log-log link.

Value

A list containing the objects MODEL, DATA, FIXED, MCMC and call. The MCMC samples are stored in the object MCMC as follows:

• z A matrix containing the MCMC samples for the spatial random field. Each column is one sample.

• mu A matrix containing the MCMC samples for the mean response (a transformation of z). Each column is one sample.

• beta A matrix containing the MCMC samples for the regressor coefficients. Each column is one sample.

• ssq A vector with the MCMC samples for the partial

• phi A vector with the MCMC samples for the spatial range parameter, if sampled.

• omg A vector with the MCMC samples for the relative nugget parameter, if sampled.

• logLik A vector containing the value of the log-likelihood evaluated at each sample.

• acc_ratio The acceptance ratio for the joint update of the parameters phi and omg, if sampled.

• sys_time The total computing time for the MCMC sampling.

• Nout, Nbi, Nthin As in input. Used internally in other functions.

The other objects contain input variables. The object call contains the function call.

Examples

## Not run: 
data(rhizoctonia)

### Create prediction grid
predgrid <- mkpredgrid2d(rhizoctonia[c("Xcoord", "Ycoord")],
                         par.x = 100, chull = TRUE, exf = 1.2)

### Combine observed and prediction locations
rhizdata <- stackdata(rhizoctonia, predgrid$grid)
##'
### Define the model
corrf <- "spherical"
family <- "binomial.probit"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
phiprior <- c(100, 1, 1000, 100) # U(100, 200)
phisc <- 3
omgprior <- c(2, 1, 1, 0)        # Exp(mean = 2)
omgsc <- .1
##'
### MCMC sizes
Nout <- 100
Nthin <- 1
Nbi <- 0

### Trial run
emt <- mcsglmm(Infected ~ 1, family, rhizdata, weights = Total,
               atsample = ~ Xcoord + Ycoord,
               Nout = Nout, Nthin = Nthin, Nbi = Nbi,
               betm0 = betm0, betQ0 = betQ0, ssqdf = ssqdf, ssqsc = ssqsc,
               corrpriors = list(phi = phiprior, omg = omgprior), 
               corrfcn = corrf, kappa = kappa,
               corrtuning = list(phi = phisc, omg = omgsc, kappa = 0),
               dispersion = 1, test = 10)

### Full run
emc <- update(emt, test = FALSE)

emcmc <- mcmcmake(emc)
summary(emcmc[, c("phi", "omg", "beta", "ssq")])
plot(emcmc[, c("phi", "omg", "beta", "ssq")])

## End(Not run)

MCMC samples from the Spatial GLMM

Description

Draw MCMC samples from the Spatial GLMM with known link function

Usage

mcsglmm_mala(
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  linkp,
  phi,
  omg,
  kappa,
  Nout,
  Nthin = 1,
  Nbi = 0,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  corrpriors,
  corrtuning,
  malatuning,
  dispersion = 1,
  longlat = FALSE,
  test = FALSE
)

Arguments

Details

The four-parameter prior for phi is defined by

\propto (\phi - \theta_4)^{\theta_2 -1} \exp\{-(\frac{\phi - \theta_4}{\theta_1})^{\theta_3}\}

for \phi > \theta_4. The prior for omg is similar. The prior parameters correspond to scale, shape, exponent, and location. See arXiv:1005.3274 for details of this distribution.

The GEV (Generalised Extreme Value) link is defined by

\mu = 1 - \exp\{-\max(0, 1 + \nu x)^{\frac{1}{\nu}}\}

for any real \nu. At \nu = 0 it reduces to the complementary log-log link.

Value

A list containing the objects MODEL, DATA, FIXED, MCMC and call. The MCMC samples are stored in the object MCMC as follows:

• z A matrix containing the MCMC samples for the spatial random field. Each column is one sample.

• mu A matrix containing the MCMC samples for the mean response (a transformation of z). Each column is one sample.

• beta A matrix containing the MCMC samples for the regressor coefficients. Each column is one sample.

• ssq A vector with the MCMC samples for the partial

• phi A vector with the MCMC samples for the spatial range parameter, if sampled.

• omg A vector with the MCMC samples for the relative nugget parameter, if sampled.

• logLik A vector containing the value of the log-likelihood evaluated at each sample.

• acc_ratio The acceptance ratio for the joint update of the parameters phi and omg, if sampled.

• sys_time The total computing time for the MCMC sampling.

• Nout, Nbi, Nthin As in input. Used internally in other functions.

The other objects contain input variables. The object call contains the function call.

Examples

## Not run: 
data(rhizoctonia)

### Create prediction grid
predgrid <- mkpredgrid2d(rhizoctonia[c("Xcoord", "Ycoord")],
                         par.x = 100, chull = TRUE, exf = 1.2)

### Combine observed and prediction locations
rhizdata <- stackdata(rhizoctonia, predgrid$grid)
##'
### Define the model
corrf <- "spherical"
family <- "binomial.probit"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
phiprior <- c(100, 1, 1000, 100) # U(100, 200)
phisc <- 3
omgprior <- c(2, 1, 1, 0)        # Exp(mean = 2)
omgsc <- .1
##'
### MCMC sizes
Nout <- 100
Nthin <- 1
Nbi <- 0

### Trial run
emt <- mcsglmm_mala(Infected ~ 1, family, rhizdata, weights = Total,
               atsample = ~ Xcoord + Ycoord,
               Nout = Nout, Nthin = Nthin, Nbi = Nbi,
               betm0 = betm0, betQ0 = betQ0, ssqdf = ssqdf, ssqsc = ssqsc,
               corrpriors = list(phi = phiprior, omg = omgprior), 
               corrfcn = corrf, kappa = kappa,
               corrtuning = list(phi = phisc, omg = omgsc, kappa = 0),
               malatuning = .003, dispersion = 1, test = 10)

### Full run
emc <- update(emt, test = FALSE)

emcmc <- mcmcmake(emc)
summary(emcmc[, c("phi", "omg", "beta", "ssq")])
plot(emcmc[, c("phi", "omg", "beta", "ssq")])

## End(Not run)

MCMC samples from the transformed Gaussian model

Description

Draw MCMC samples from the transformed Gaussian model with known link function

Usage

mcstrga(
  formula,
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  linkp,
  phi,
  omg,
  kappa,
  Nout,
  Nthin = 1,
  Nbi = 0,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  tsqdf,
  tsqsc,
  corrpriors,
  corrtuning,
  longlat = FALSE,
  test = FALSE
)

Arguments

Details

Simulates from the posterior distribution of this model.

Value

A list containing the objects MODEL, DATA, FIXED, MCMC and call. The MCMC samples are stored in the object MCMC as follows:

• z A matrix containing the MCMC samples for the spatial random field. Each column is one sample.

• mu A matrix containing the MCMC samples for the mean response (a transformation of z). Each column is one sample.

• beta A matrix containing the MCMC samples for the regressor coefficients. Each column is one sample.

• ssq A vector with the MCMC samples for the partial

• tsq A vector with the MCMC samples for the measurement error variance.

• phi A vector with the MCMC samples for the spatial range parameter, if sampled.

• omg A vector with the MCMC samples for the relative nugget parameter, if sampled.

• logLik A vector containing the value of the log-likelihood evaluated at each sample.

• acc_ratio The acceptance ratio for the joint update of the parameters phi and omg, if sampled.

• sys_time The total computing time for the MCMC sampling.

• Nout, Nbi, Nthin As in input. Used internally in other functions.

The other objects contain input variables. The object call contains the function call.

Examples

## Not run: 
### Load the data
data(rhizoctonia)
rhiz <- na.omit(rhizoctonia)
rhiz$IR <- rhiz$Infected/rhiz$Total # Incidence rate of the
                              # rhizoctonia disease

### Define the model
corrf <- "spherical"
ssqdf <- 1
ssqsc <- 1
tsqdf <- 1
tsqsc <- 1
betm0 <- 0
betQ0 <- diag(.01, 2, 2)
phiprior <- c(200, 1, 1000, 100) # U(100, 300)
phisc <- 1
omgprior <- c(3, 1, 1000, 0) # U(0, 3)
omgsc <- 1
linkp <- 1

## MCMC parameters
Nout <- 100
Nbi <- 0
Nthin <- 1

samplt <- mcstrga(Yield ~ IR, data = rhiz,
                  atsample = ~ Xcoord + Ycoord, corrf = corrf,
                  Nout = Nout, Nthin = Nthin,
                  Nbi = Nbi, betm0 = betm0, betQ0 = betQ0,
                  ssqdf = ssqdf, ssqsc = ssqsc,
                  tsqdf = tsqdf, tsqsc = tsqsc,
                  corrprior = list(phi = phiprior, omg = omgprior),
                  linkp = linkp,
                  corrtuning = list(phi = phisc, omg = omgsc, kappa = 0),
                  test=10)

sample <- update(samplt, test = FALSE)

## End(Not run)

MCMC samples from the transformed Gaussian model

Description

Draw MCMC samples from the transformed Gaussian model with known link function

Usage

mcstrga_mala(
  formula,
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  linkp,
  phi,
  omg,
  kappa,
  Nout,
  Nthin = 1,
  Nbi = 0,
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  tsqdf,
  tsqsc,
  corrpriors,
  corrtuning,
  malatuning,
  longlat = FALSE,
  test = FALSE
)

Arguments

Details

Simulates from the posterior distribution of this model.

Value

A list containing the objects MODEL, DATA, FIXED, MCMC and call. The MCMC samples are stored in the object MCMC as follows:

• z A matrix containing the MCMC samples for the spatial random field. Each column is one sample.

• mu A matrix containing the MCMC samples for the mean response (a transformation of z). Each column is one sample.

• beta A matrix containing the MCMC samples for the regressor coefficients. Each column is one sample.

• ssq A vector with the MCMC samples for the partial

• tsq A vector with the MCMC samples for the measurement error variance.

• phi A vector with the MCMC samples for the spatial range parameter, if sampled.

• omg A vector with the MCMC samples for the relative nugget parameter, if sampled.

• logLik A vector containing the value of the log-likelihood evaluated at each sample.

• acc_ratio The acceptance ratio for the joint update of the parameters phi and omg, if sampled.

• sys_time The total computing time for the MCMC sampling.

• Nout, Nbi, Nthin As in input. Used internally in other functions.

The other objects contain input variables. The object call contains the function call.

Examples

## Not run: 
### Load the data
data(rhizoctonia)
rhiz <- na.omit(rhizoctonia)
rhiz$IR <- rhiz$Infected/rhiz$Total # Incidence rate of the
                              # rhizoctonia disease

### Define the model
corrf <- "spherical"
ssqdf <- 1
ssqsc <- 1
tsqdf <- 1
tsqsc <- 1
betm0 <- 0
betQ0 <- diag(.01, 2, 2)
phiprior <- c(200, 1, 1000, 100) # U(100, 300)
phisc <- 1
omgprior <- c(3, 1, 1000, 0) # U(0, 3)
omgsc <- 1
linkp <- 1

## MCMC parameters
Nout <- 100
Nbi <- 0
Nthin <- 1

samplt <- mcstrga_mala(Yield ~ IR, data = rhiz,
                  atsample = ~ Xcoord + Ycoord, corrf = corrf,
                  Nout = Nout, Nthin = Nthin,
                  Nbi = Nbi, betm0 = betm0, betQ0 = betQ0,
                  ssqdf = ssqdf, ssqsc = ssqsc,
                  tsqdf = tsqdf, tsqsc = tsqsc,
                  corrprior = list(phi = phiprior, omg = omgprior),
                  linkp = linkp,
                  corrtuning = list(phi = phisc, omg = omgsc, kappa = 0),
                  malatuning = .0002, test=10)

sample <- update(samplt, test = FALSE)

## End(Not run)

Make prediction grid

Description

This function creates a grid for prediction.

Usage

mkpredgrid2d(
  pnts.x,
  pnts.y,
  par.x,
  par.y,
  isby = FALSE,
  chull = FALSE,
  exf = 1
)

Arguments

Details

If chull this function first calculates the convex hull of the points. If exf is not 1 the borders are expanded. Then the function calls point.in.polygon to select points that fall inside the borders.

Value

A list with components

• grid A two-column matrix with the prediction grid.

• xycoord A list with components "x" and "y" containing the sequence of points used to create the grid.

• xygrid A matrix with the full square grid derived from xycoord.

• borders The (expanded) borders of the domain.

• inxygrid A logical vector indicating which rows of xycoord fall inside borders, and therefore correspond to the grid.

Examples

## Not run: 
data(rhizoctonia)
predgrid <- mkpredgrid2d(rhizoctonia[c("Xcoord", "Ycoord")],
                         par.x = 100, chull = TRUE, exf = 1.2)
plot(predgrid$borders, type = "l")         # Domain for prediction
points(predgrid$grid, pch = 20, cex = .3)  # Prediction locations
points(rhizoctonia[c("Xcoord", "Ycoord")]) # Observed locations

## End(Not run)

Plot the estimated Bayes factors

Description

This function plots the estimated logarithm Bayes factors from the function bf2new.

Usage

plotbf2(
  bf2obj,
  pars = c("linkp", "phi", "omg", "kappa"),
  profile = length(pars) > 2,
  ...
)

Arguments

Details

Depending on whether pars has length 1 or 2, this function creates a line or a contour plot of the estimated Bayes factors. If its length is 3 or 4, then it produces multiple profile plots. In this case the variable is fixed at different values and the maximum Bayes factor corresponding to the fixed value is plotted against that value.

Value

This function returns nothing.

Examples

## Not run: 
data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
family <- "binomial.probit"
### Skeleton points
philist <- c(100, 140, 180)
omglist <- c(.5, 1)
parlist <- expand.grid(linkp=0, phi=philist, 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, family, rhizoctonia, weights = Total,
                       atsample = ~ Xcoord + Ycoord,
                       Nout = Nout, Nthin = Nthin, Nbi = Nbi,
                       betm0 = betm0, betQ0 = betQ0,
                       ssqdf = ssqdf, ssqsc = ssqsc,
                       phi = parlist$phi[i], omg = parlist$omg[i],
                       linkp = parlist$linkp[i], kappa = parlist$kappa[i],
                       corrfcn = corrf,
                       corrtuning=list(phi = 0, omg = 0, kappa = 0))
}
bf <- bf1skel(runs)
bfall <- bf2new(bf, phi = seq(100, 200, 10), omg = seq(0, 2, .2))
plotbf2(bfall, c("phi", "omg"))
plotbf2(bfall, c("phi", "omg"), profile = TRUE, type = "b", ylab="log(BF)")

## End(Not run)

Reverse logistic regression estimation

Description

Perform the reverse logistic regression estimation

Usage

revlogreg(lglk, N)

Arguments

Details

Estimation is done by maximising the reverse logistic log likelihood.

Value

A vector containing the reverse logistic regression estimates of the logarithm of the Bayes factors. The first set of parameters is taken as the reference model so its estimate is always 0.

References

Geyer, C. J. (1994). Estimating normalizing constants and reweighting mixtures in Markov chain Monte Carlo. Technical Report 568, School of Statistics, University of Minnesota.

Rhizoctonia root rot infections

Description

Rhizoctonia root rot infections.

Usage

data(rhizoctonia)

Format

A data frame with 100 rows and 5 variables.

Details

A dataset containing the number of infected roots and the sample coordinate. The data were collected by Dr Jim Cook at Washington State University.

• Xcoord Longitude of the sampling location.

• Ycoord Latitude of the sampling location.

• Total Total number of roots sampled at that location.

• Infected Number of infected roots found at that location.

• Yield Barley yield at that location. These data were obtained by hand-harvesting a 4-square-meter area in the sampling location. One observation is missing.

Note

We acknowledge Hao Zhang for providing these data.

Source

See reference below.

References

Zhang, H. (2002). On estimation and prediction for spatial generalized linear mixed models. Biometrics, 58(1), 129-136.

Simulation from a spatial model

Description

Simulate from a variety of spatial models.

Usage

rsglmm(
  n,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  beta,
  linkp,
  phi,
  omg,
  kappa,
  ssq,
  corrfcn = "matern",
  longlat = FALSE,
  dispersion = 1,
  returnGRF = FALSE,
  warndisp = TRUE
)

rstrga(
  n,
  formula,
  data,
  weights,
  subset,
  offset,
  atsample,
  beta,
  linkp,
  phi,
  omg,
  kappa,
  ssq,
  corrfcn = "matern",
  longlat = FALSE,
  dispersion = 1,
  returnGRF = FALSE
)

rsgrf(
  n,
  formula,
  data,
  subset,
  offset,
  atsample,
  beta,
  phi,
  omg,
  kappa,
  ssq,
  corrfcn = "matern",
  longlat = FALSE
)

Arguments

Details

The spatial Gaussian random field is simulated using the Cholesky decomposition of the covariance matrix.

The sample from a quasi distribution uses a hack which matches the mean and the variance of the distribution. See the source code for details.

Value

A data frame containing the predictors, sampling locations, optional weights, and samples.

Examples

## Not run: 
n <- 100
beta <- c(-2, 1)
phi <- .2
omg <- .3
linkp <- 0
ssq <- 1
l <- rep(10, n)
corrf <- "matern"
kappa <- .5
family <- "poisson"
Xcoord <- runif(n)
Ycoord <- runif(n)
f <- Xcoord + Ycoord
formula <- y|z ~ f
mydata <- rsglmm(1, formula, family, weights = l,
                 atsample = ~ Xcoord + Ycoord, beta = beta, linkp = linkp,
                 phi = phi, omg = omg, kappa = kappa, ssq = ssq,
                 corrfcn = corrf, returnGRF = TRUE)

## End(Not run)

Selection of multiple importance sampling distributions

Description

Selection of multiple importance sampling distributions

Usage

select_proposals_SEQ(
  pargrid,
  K,
  istart,
  relativeSE = FALSE,
  N1,
  N2,
  Nthin,
  Nbi,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  dispersion = 1,
  longlat = FALSE,
  nbatch1 = 0.5,
  nbatch2 = 0.5,
  S1method = c("RL", "MW"),
  bvmethod = c("Standard", "TukeyHanning", "Bartlett"),
  transf = c("no", "mu", "wo")
)

select_proposals_MNX(
  pargrid,
  istart,
  nfix,
  relativeSE = FALSE,
  N1,
  N2,
  Nthin,
  Nbi,
  cooling,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  dispersion = 1,
  longlat = FALSE,
  nbatch1 = 0.5,
  nbatch2 = 0.5,
  S1method = c("RL", "MW"),
  bvmethod = c("Standard", "TukeyHanning", "Bartlett"),
  transf = c("no", "mu", "wo"),
  verbose = FALSE
)

select_proposals_ENT(
  pargrid,
  istart,
  nfix,
  relativeSE = FALSE,
  N1,
  Nthin,
  Nbi,
  cooling,
  formula,
  family = "gaussian",
  data,
  weights,
  subset,
  offset,
  atsample,
  corrfcn = "matern",
  betm0,
  betQ0,
  ssqdf,
  ssqsc,
  dispersion = 1,
  longlat = FALSE,
  nbatch1 = 0.5,
  nbatch2 = 0.5,
  S1method = c("RL", "MW"),
  bvmethod = c("Standard", "TukeyHanning", "Bartlett"),
  transf = c("no", "mu", "wo"),
  verbose = FALSE
)

Arguments

Details

SEQ

is a sequential method starting with istart and additng to it until K proposals have been selected. At each iteration, the point with the highest (relative?) standard error is added

MNX

is the minimax method. The chosen proposal corresponds to the lowest maximum (relative?) standard error.

ENT

is the entropy method. The chosen proposal corresponds to the highest determinant of the (relative?) covariance matrix at the first stage.

Value

A list with components

selected

The rows of pargrid selected.

isel

The indices of the rows of pargrid selected.

se

The standard error corresponding to the selected parameters.

samples

A list containing the samples from the selected parameters.

References

Roy, V., & Evangelou, E. (2018). Selection of proposal distributions for generalized importance sampling estimators. arXiv preprint arXiv:1805.00829.

Examples

## Not run: 
data(rhizoctonia)
### Define the model
corrf <- "spherical"
kappa <- 0
ssqdf <- 1
ssqsc <- 1
betm0 <- 0
betQ0 <- .01
family <- "binomial.probit"
formula <- Infected ~ 1
atsample <- ~ Xcoord + Ycoord
### Skeleton points
philist <- seq(100, 200, 10)
omglist <- 0
parlist <- expand.grid(linkp=0, phi=philist, omg=omglist, kappa = kappa)
### MCMC sizes
Nout <- 100
Nthin <- 1
Nbi <- 10
## Select proposals
K <- 3                        # Choose 3 proposals
istart_SEQ <- 6               # Start with middle
istart_MNX <- istart_ENT <- c(6, 2, 10)
cooling_MNX <- .05/log((0:24 %/% 5)*5 + exp(1))
cooling_ENT <- .3/log((0:49 %/% 10)*10 + exp(1))
prop_SEQ <- select_proposals_SEQ(pargrid = parlist, K = K,
                                 istart = istart_SEQ,
                                 relativeSE = TRUE, 
                                 N1 = Nout, N2 = Nout,
                                 Nthin = Nthin, Nbi = Nbi,
                                 formula = formula, family = family,
                                 data = rhizoctonia, weights = Total, 
                                 atsample = atsample, corrfcn = corrf,
                                 betm0 = betm0, betQ0 = betQ0,
                                 ssqdf = ssqdf, ssqsc = ssqsc,
                                 dispersion = 1, longlat = FALSE,
                                 nbatch1 = 0.5, nbatch2 = 0.5,
                                 bvmethod = "TukeyHanning",
                                 transf = "mu")
prop_MNX <- select_proposals_MNX(pargrid = parlist,
                                 istart = istart_MNX, nfix = 1L,
                                 cooling = cooling_MNX, 
                                 relativeSE = TRUE, 
                                 N1 = Nout, N2 = Nout,
                                 Nthin = Nthin, Nbi = Nbi,
                                 formula = formula, family = family,
                                 data = rhizoctonia, weights = Total, 
                                 atsample = atsample, corrfcn = corrf,
                                 betm0 = betm0, betQ0 = betQ0,
                                 ssqdf = ssqdf, ssqsc = ssqsc,
                                 dispersion = 1, longlat = FALSE,
                                 nbatch1 = 0.5, nbatch2 = 0.5,
                                 bvmethod = "TukeyHanning",
                                 transf = "mu",
                                 verbose = TRUE)
prop_ENT <- select_proposals_ENT(pargrid = parlist,
                                 istart = istart_ENT, nfix = 1L,
                                 cooling = cooling_ENT, 
                                 relativeSE = TRUE, 
                                 N1 = Nout, 
                                 Nthin = Nthin, Nbi = Nbi,
                                 formula = formula, family = family,
                                 data = rhizoctonia, weights = Total, 
                                 atsample = atsample, corrfcn = corrf,
                                 betm0 = betm0, betQ0 = betQ0,
                                 ssqdf = ssqdf, ssqsc = ssqsc,
                                 dispersion = 1, longlat = FALSE,
                                 nbatch1 = 0.5, nbatch2 = 0.5,
                                 bvmethod = "TukeyHanning",
                                 transf = "mu",
                                 verbose = TRUE)

## End(Not run)

Spatial variance-covariance matrix

Description

Calculates the spatial variance-covariance matrix for a selection of correlation functions.

Usage

spcovariance(...)

## S3 method for class 'formula'
spcovariance(
  formula,
  data,
  subset,
  corrfcn,
  ssq,
  phi,
  omg,
  kappa,
  longlat = FALSE,
  ...
)

## S3 method for class 'numeric'
spcovariance(x, corrfcn, ssq, phi, omg, kappa, ...)

## S3 method for class 'dist'
spcovariance(x, corrfcn, ssq, phi, omg, kappa, ...)

Arguments

Value

For a formula input, a variance-covariance matrix. For a numeric input, an object of the the same dimensions as its first input.

Spatial log likelihood

Description

Spatial joint log likelihood

Usage

sploglik(pargrid, runs, transf = c("no", "mu", "wo"))

Arguments

Details

Computes the joint log likelihood log f(y,T(z)|parameters) where T(z) is the transformation, for each (y,z) in runs and for parameters in pargrid up to a constant which does not depend on the parameters. The parameters beta and sigma^2 are integrated out.

Value

A matrix with number of rows the total number of samples in runs and number of columns the number of rows in pargrid. The [i,j] element of the matrix is the value of the loglikelihood at the ith sample when all samples in runs are put together evaluated at the jth parameter value.

Spatial log likelihood

Description

Spatial joint log likelihood

Usage

sploglik_cross(runs, transf = c("no", "mu", "wo"))

Arguments

Details

Computes the joint log likelihood log f(y,T(z)|parameters) where T(z) is the transformation, for each (y,z) in runs and for parameters in runs up to a constant which does not depend on the parameters. The parameters beta and sigma^2 are integrated out.

Value

A matrix with number of rows the total number of samples in runs and number of columns the length of runs. The [i,j] element of the matrix is the value of the loglikelihood at the ith sample when all samples in runs are put together evaluated at the jth parameter value.

Combine data.frames

Description

Combine data.frames

Usage

stackdata(..., fillwith = NA, keepclass = FALSE)

Arguments

Details

This function combines data.frames by filling in missing variables. This is useful for combining data from sampled locations with prediction locations.

If fillwith is a named object, its names must correspond to the names of variables in the data frames. If a variable is missing, then it is filled with the corresponding value in fillwith. fillwith can contain only one unnamed component which corresponds to the default filling.

Value

A stacked data.frame.

Examples

## Not run: 
d1 <- data.frame(w = 1:3, z = 4:6 + 0.1)
d2 <- data.frame(w = 3:7, x = 1:5, y = 6:10)
(d12a <- stackdata(d1, d2))
lapply(d12a, class)
(d12b <- stackdata(d1, d2, fillwith = c(x = NA, y = 0, z = -99)))
lapply(d12b, class)
(d12c <- stackdata(d1, d2, fillwith = c(x = NA, y = 0, z = -99),
                   keepclass = TRUE))
lapply(d12c, class)
(d12d <- stackdata(d1, d2, fillwith = c(x = NA, 0)))

data(rhizoctonia)
predgrid <- mkpredgrid2d(rhizoctonia[c("Xcoord", "Ycoord")],
                         par.x = 100, chull = TRUE, exf = 1.2)
rhizdata <- stackdata(rhizoctonia, predgrid$grid)

## End(Not run)

Subset MCMC chain

Description

Return subset of MCMC chain.

Usage

## S3 method for class 'geomcmc'
subset(x, subset, ...)

Arguments

Value

A similar object as x with the subsetted chain.