# MCMC samples from the Spatial GLMM Draw MCMC samples from the Spatial GLMM with known link function ```r 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 - `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. - `offset`: See `lm`. - `atsample`: A formula in the form `~ x1 + x2 + ... + xd` with the coordinates of the sampled locations. - `corrfcn`: Spatial correlation function. See `geoBayes_correlation` for details. - `linkp`: Parameter of the link function. A scalar value. - `phi`: Optional starting value for the MCMC for the spatial range parameter `phi`. Defaults to the mean of its prior. If `corrtuning[["phi"]]` is 0, then this argument is required and it corresponds to the fixed value of `phi`. This can be a vector of the same length as Nout. - `omg`: Optional starting value for the MCMC for the relative nugget parameter `omg`. Defaults to the mean of its prior. If `corrtuning[["omg"]]` is 0, then this argument is required and it corresponds to the fixed value of `omg`. This can be a vector of the same length as Nout. - `kappa`: Optional starting value for the MCMC for the spatial correlation parameter `kappa` (Matern smoothness or exponential power). Defaults to the mean of its prior. If `corrtuning[["kappa"]]` is 0 and it is needed for the chosen correlation function, then this argument is required and it corresponds to the fixed value of `kappa`. This can be a vector of the same length as Nout. - `Nout`: Number of MCMC samples to return. This can be a vector for running independent chains. 0 elements are dropped. - `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. - `corrpriors`: A list with the components `phi`, `omg` and `kappa` as needed. These correspond to the prior distribution parameters. For `phi` and `omg` it must be a vector of length 4. The generalized inverse gamma prior is assumed and the input corresponds to the parameters scale, shape, exponent, location in that order (see Details). For `kappa` it must be a vector of length 2. A uniform prior is assumed and the input corresponds to the lower and upper bounds in that order. - `corrtuning`: A vector or list with the components `phi`, `omg` and `kappa` as needed. These correspond to the random walk parameter for the Metropolis-Hastings step. Smaller values increase the acceptance ratio. Set this to 0 for fixed parameter value. - `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. ## Returns 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. ## 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}\}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)^{\frac{1}{\nu}}\}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. ## Examples ```r ## 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) ```

mcsglmm function