# mcsglmm: MCMC samples from the Spatial GLMM In geoBayes: Analysis of Geostatistical Data using Bayes and Empirical Bayes Methods

## Description

Draw MCMC samples from the Spatial GLMM with known link function

## Usage

 ```1 2 3 4 5 6``` ```mcsglmm(formula, family = c("gaussian", "binomial", "poisson", "Gamma", "GEV.binomial", "GEVD.binomial", "Wallace.binomial"), data, weights, subset, atsample, Nout, Nthin = 1, Nbi = 0, betm0, betQ0, ssqdf, ssqsc, phipars, omgpars, corrfcn = c("matern", "spherical", "powerexponential"), kappa, linkp, phisc, omgsc, zstart, phistart, omgstart, dispersion = 1, longlat = FALSE, test = FALSE) ```

## Arguments

 `formula` A representation of the model in the form `response ~ terms`. The response must be set to `NA`'s at the prediction locations (see the examples on how to do this using the function `stackdata`). At the observed locations the response is assumed to be a total of replicated measurements. The number of replications is inputted using the argument `weights`. `family` The distribution of the data. The `"GEVbinomial"` family is the binomial family with link the GEV link (see Details). `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 `~ x1 + x2 + ... + xd` with the coordinates of the sampled locations. `Nout` Number of MCMC samples to return. `Nthin` The thinning of the MCMC algorithm. `Nbi` The burn-in of the MCMC algorithm. `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 chi-square prior for the partial sill parameter. `ssqsc` Scale for the scaled inverse chi-square prior for the partial sill parameter. `phipars` Parameters for the generalized inverse gamma prior for the spatial range parameter `phi`. A four dimensional vector with parameters scale, shape, exponent, location in that order. See Details. `omgpars` Parameters for the generalized inverse gamma prior for the relative nugget parameter `omg`. A four dimensional vector with parameters scale, shape, exponent, location in that order. See Details. `corrfcn` Spatial correlation function. See `ebsglmm` for details. `kappa` Spatial correlation parameter. Smoothness parameter for Matern, exponent for the power family. `linkp` Parameter of the link function. For binomial, a positive number for the degrees of freedom of the robit family or "logit" or "probit". For the other families any number for the exponent of the Box-Cox transformation. `phisc` Random walk parameter for `phi`. Smaller values increase the acceptance ratio. Set this to 0 for fixed `phi`. In this case the fixed value is given in the argument `phistart`. `omgsc` Random walk parameter for `omg`. Smaller values increase the acceptance ratio. Set this to 0 for fixed `omg`. In this case the fixed value is given in the argument `omgstart`. `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. `phistart` Optional starting value for the MCMC for the spatial range parameter `phi`. Defaults to the mean of its prior. If `phisc` is 0, then this argument is required and it corresponds to the fixed value of `phi`. `omgstart` Optional starting value for the MCMC for the relative nugget parameter `omg`. Defaults to the mean of its prior. If `omgsc` is 0, then this argument is required and itcorresponds to the fixed value of `omg`. `dispersion` The fixed dispersion parameter. `longlat` How to compute the distance between locations. If `FALSE`, Euclidean distance, if `TRUE` Great Circle distance. See `spDists`. `test` Whether this is a trial run to monitor the acceptance ratio of the random walk for `phi` and `omg`. If set to `TRUE`, the acceptance ratio will be printed on the screen every 100 iterations of the MCMC. Tune the `phisc` and `omgsc` parameters in order to achive 20 to 30% acceptance. Set this to a positive number to change the default 100. No thinning or burn-in are done when testing.

## Details

The four-parameter prior for `phi` is defined by

propto (phi - phiprior[4])^(phiprior[2]-1) * exp(-((phi-phiprior[4])/phiprior[1])^phiprior[3])

for phi > phiprior[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)^(1/nu)]

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

## Value

A list containing the MCMC samples and other variables 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 sill parameter.

• `phi` A vector with the MCMC samples for the spatial range parameter.

• `omg` A vector with the MCMC samples for the relative nugget parameter.

• `nu` The link function parameter translated to numeric code used internally.

• `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`.

• `sys_time` The total computing time for the MCMC sampling.

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

• `response` The value of the response variable at the observed locations. Used internally in other functions.

• `weights` The response weights at the observed locations. Used internally in other functions.

• `modelmatrix` The model matrix at the observed locations. Used internally in other functions.

• `family` As in input. Used internally in other functions.

• `betm0`, `betQ0`, `ssqdf`, `ssqsc`, `corrfcn`, `kappa`, `dispersion` As in input. Used internally in other functions.

• `locations` Coordinates of the observed locations. Used internally in other functions.

• `whichobs` A logical vector indicated which rows in the data and in the MCMC samples for the spatial random field correspond to the observed locations.

## 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 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54``` ```## 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" 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 linkp <- "probit" ### MCMC sizes Nout <- 100 Nthin <- 1 Nbi <- 0 ### Trial run emt <- mcsglmm(Infected ~ 1, 'binomial', rhizdata, weights = Total, atsample = ~ Xcoord + Ycoord, Nout = Nout, Nthin = Nthin, Nbi = Nbi, betm0 = betm0, betQ0 = betQ0, ssqdf = ssqdf, ssqsc = ssqsc, phipars = phiprior, omgpars = omgprior, linkp = linkp, corrfcn = corrf, kappa = kappa, phisc = phisc, omgsc = omgsc, dispersion = 1, test = 10) ### Full run emc <- mcsglmm(Infected ~ 1, 'binomial', rhizdata, weights = Total, atsample = ~ Xcoord + Ycoord, Nout = Nout, Nthin = Nthin, Nbi = Nbi, betm0 = betm0, betQ0 = betQ0, ssqdf = ssqdf, ssqsc = ssqsc, phipars = phiprior, omgpars = omgprior, linkp = linkp, corrfcn = corrf, kappa = kappa, phisc = phisc, omgsc = omgsc, dispersion = 1, test = FALSE) plot.ts(cbind(phi = emc\$phi, omg = emc\$omg, beta = c(emc\$beta), ssq = emc\$ssq), nc = 2) emcmc <- mcmcmake(emc) summary(emcmc[, c("phi", "omg", "beta", "ssq")]) ## End(Not run) ```

geoBayes documentation built on May 30, 2017, 5:47 a.m.