Empirical Bayes estimation for the TGRFM
Description
Empirical Bayes estimation for the spatial transformed Gaussian model.
Usage
1 2 3 4 5 6  ebstrga(formula, data, weights, subset, atsample, parskel, paroptim,
corrfcn = c("matern", "spherical", "powerexponential"), Nout, Nthin = 1,
Nbi = 0, Npro, Nprt = 1, Nprb = 0, betm0, betQ0, ssqdf, ssqsc, tsqdf,
tsqsc, zstart, dispersion = 1, bfsize1 = 0.8, reference = 1,
bfmethod = c("RL", "MW"), transf = FALSE, useCV = TRUE,
longlat = FALSE, control = list(), verbose = TRUE)

Arguments
formula 
A representation of the model in the form

data 
An optional data frame containing the variables in the model. 
weights 
An optional vector of weights. Number of replicated samples for Gaussian and gamma, number of trials for binomial, time length for Poisson. 
subset 
An optional vector specifying a subset of observations to be used in the fitting process. 
atsample 
A formula in the form 
parskel 
A data frame with the components "linkp", "phi",
"omg", and "kappa", corresponding to the link function, the
spatial range, the relative nugget, and the spatial smoothness
parameters. The latter can be omitted if not used in the
correlation function. Let k denote the number of rows. Then, k
different MCMC samples will be taken from the models with
parameters fixed at those values. For a square grid the output
from the function 
paroptim 
A named list with the components "linkp", "phi", "omg", "kappa". The latter can be omitted if not used in the correlation function. Each component must be numeric with length 1, 2, or 3 with elements in increasing order but for the binomial family linkp is also allowed to be the character "logit" and "probit". If its length is 1, then the corresponding parameter is considered to be fixed at that value. If 2, then the two numbers denote the lower and upper bounds for the optimisation of that parameter (infinities are allowed). If 3, these correspond to lower bound, starting value, upper bound for the estimation of that parameter. 
corrfcn 
Spatial correlation function. See

Nout 
A scalar or vector of size k. Number of MCMC samples
to take for each run of the MCMC algorithm for the estimation of
the Bayes factors. See argument 
Nthin 
A scalar or vector of size k. The thinning of the MCMC algorithm for the estimation of the Bayes factors. 
Nbi 
A scalar or vector of size k. The burnin of the MCMC algorithm for the estimation of the Bayes factors. 
Npro 
A scalar. The number of Gibbs samples to take for estimation of the conjugate parameters and for prediction at the unsampled locations while the other parameters are fixed at their empirical Bayes estimates. 
Nprt 
The thinning of the Gibbs algorithm for the estimation of the conjugate parameters and for prediction. 
Nprb 
The burnin of the Gibbs algorithm for the estimation of the conjugate parameters and for prediction. 
betm0 
Prior mean for beta (a vector or scalar). 
betQ0 
Prior standardised precision (inverse variance) matrix. Can be a scalar, vector or matrix. The first two imply a diagonal with those elements. Set this to 0 to indicate a flat improper prior. 
ssqdf 
Degrees of freedom for the scaled inverse chisquare prior for the partial sill parameter. 
ssqsc 
Scale for the scaled inverse chisquare prior for the partial sill parameter. 
tsqdf 
Degrees of freedom for the scaled inverse chisquare prior for the measurement error parameter. 
tsqsc 
Scale for the scaled inverse chisquare prior for the measurement error parameter. 
zstart 
Optional starting value for the MCMC for the GRF.
This can be either a scalar, a vector of size n where n is the
number of sampled locations, or a matrix with dimensions n by k
where k is the number of the skeleton points in 
dispersion 
The fixed dispersion parameter. 
bfsize1 
A scalar or vector of the same length as 
reference 
Which model goes in the denominator of the Bayes factors. 
bfmethod 
Which method to use to calculate the Bayes factors: Reverse logistic or MengWong. 
transf 
Whether to use the transformed sample mu for the computations. Otherwise it uses z. 
useCV 
Whether to use control variates for finer corrections. 
longlat 
How to compute the distance between locations. If

control 
A list of control parameters for the optimisation.
See 
verbose 
Whether to print messages when completing each stage on screen. 
Details
Runs the MCMC sampling, computes the importance weights, and estimates the parameters.
Value
A list with components

parest
The parameter paroptims 
skeleton
The skeleton points used with the corresponding logarithm of the Bayes factors at those points. 
optim
The output from theoptim
function. 
mcmcsample
The MCMC samples for the remaining parameters and the random field. These samples correspond to the Gibbs and MetropolisHasting samples after fixing the parameters estimated by empirical Bayes at their empirical Bayes estimates. 
sys_time
The time taken to complete the MCMC sampling, calculation of the importance weights, the optimization and the final MCMC sampling.
Note
To avoid numerical overflow it is recommended to scale the response variable by its geometric mean. This is strongly recommended when the exponent of the transformation is 0 or less. Numerical overflow may occur when samples which are simulated from a smaller exponent are being used to compute the likelihood corresponding to a larger exponent.
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.
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 36  ## 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)
### Skeleton points
philist < seq(120, 280, 40)
linkp < 1
omglist < seq(0, 2.5, .5)
parlist < expand.grid(phi = philist, linkp = linkp, omg = omglist)
paroptim < list(linkp = linkp, phi = c(100, 300), omg = c(0, 2))
### MCMC sizes
Nout < Npro < 100
Nthin < Nprt < 1
Nbi < Nprb < 0
### Estimate
est < ebstrga(Yield ~ IR, rhiz,
atsample = ~ Xcoord + Ycoord, parskel = parlist,
paroptim = paroptim, corrfcn = corrf,
Nout = Nout, Nthin = Nthin, Nbi = Nbi,
Npro = Npro, Nprt = Nprt, Nprb = Nprb,
betm0 = betm0, betQ0 = betQ0, ssqdf = ssqdf, ssqsc = ssqsc,
tsqdf = tsqdf, tsqsc = tsqsc,
useCV=TRUE)
est$parest
## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker. Vote for new features on Trello.