mcsglmm: MCMC samples from the Spatial GLMM

Description Usage Arguments Details Value Examples

View source: R/mcsp2.R

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:

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.

Related to mcsglmm in geoBayes...