Learn R Programming

CARBayesST (version 1.0)

ST.knorrheld.main: Fit the separable space-time Poisson log-linear model proposed by Knorr-Held and Besag (1998).

Description

The function fits the Bayesian spatio-temporal model proposed by Knorr-Held and Besag (1998) to Poisson count data. The natural log of the linear predictor is made up of 5 terms, an intercept term, two spatial random effects, correlated and independent, and two temporal random effects, correlated and independent. The independent random effects are assigned zero-mean Gaussian priors with a common variance, while the correlated effects are modelled by intrinsic conditional autoregressive models respectively. Inference is based on Markov chain Monte Carlo (McMC) simulation, using a combination of Gibbs sampling and Metropolis steps.

Usage

ST.knorrheld.main(formula, data=NULL, W, burnin=0, n.sample=1000, thin=1,  
blocksize.beta=5, prior.mean.beta=NULL, prior.var.beta=NULL, prior.tau2=NULL, 
verbose=TRUE)

Arguments

formula
A formula for the covariate part of the model, using the same notation as for the lm() function. The offsets should also be included here using the offset() function. The response and each covariate should be vectors of length (KN)*1, where K is the numbe
data
The matrix of offsets on the fitted value scale, such as the population size or the expected number of counts. Each row corresponds to an areal unit and each column to a time period.
W
A binary K by K neighbourhood matrix (where K is the number of spatial units). The jkth element equals one if areas (j, k) are spatially close (e.g. share a common border) and is zero otherwise.
burnin
The number of MCMC samples to discard as the burnin period. Defaults to 0.
n.sample
The number of MCMC samples to generate. Defaults to 1,000.
thin
The level of thinning to apply to the MCMC samples to reduce their temporal autocorrelation. Defaults to 1.
blocksize.beta
The size of the blocks in which to update the regression parameters in the MCMC algorithm. Defaults to 5.
prior.mean.beta
A vector of prior means for the regression parameters beta (Gaussian priors are assumed). Defaults to a vector of zeros.
prior.var.beta
A vector of prior variances for the regression parameters beta (Gaussian priors are assumed). Defaults to a vector with values 1000.
prior.tau2
The prior shape and scale in the form of c(shape, scale) for an Inverse-Gamma(shape, scale) prior for each random effect variance term tau2. Defaults to c(0.001, 0.001).
verbose
Logical, should the function update the user on its progress.

Value

  • formulaThe formula for the covariate and offset part of the model.
  • samplesA list containing the MCMC samples from the model.
  • fitted.valuesA vector containing the fitted value for each area and time point. The vector is ordered so that all spatial units at time one come first and then at time two and so on.
  • residualsA vector containing the residuals for each area and time point. The vector is ordered so that all spatial units at time one come first and then at time two and so on.
  • posterior.ZThis argument is NULL and is included for comparibility in the output object with the other models in this package.
  • median.ZThis argument is NULL and is included for comparibility in the output object with the other models in this package.
  • modelfitModel fit criteria including the Deviance Information Criterion (DIC), the effective number of parameters in the model (p.d), DIC3, Watanabe-Akaike Information Criterion (WAIC), and the log Marginal Predictive Likelihood (LMPL).
  • summary.resultsA table summarising some of the parameters in the model.
  • modelA text string describing the model.
  • acceptThe acceptance probabilities for the parameters.

Details

For further details about how to apply the function see the examples below and in the vignette.

References

Knorr-Held, L. and J. Besag (1998). Modelling Risk from a Disease in Time and Space. Statistics in Medicine 17, 2045-2060.

Examples

Run this code
#### Artificial data generated on a square

#### Set up a square lattice region
x.easting <- 1:10
x.northing <- 1:10
Grid <- expand.grid(x.easting, x.northing)
n <- nrow(Grid)
t <- 10


#### set up distance and neighbourhood (W, based on sharing a common border) matrices
distance <-array(0, c(n,n))
W <-array(0, c(n,n))
     for(i in 1:n)
     {
     	for(j in 1:n)
		{
		temp <- (Grid[i,1] - Grid[j,1])^2 + (Grid[i,2] - Grid[j,2])^2
		distance[i,j] <- sqrt(temp)
			if(temp==1)  W[i,j] <- 1 
		}	
	}
	
	
#### Generate data
n.all <- n * t
E <- rep(100, n.all)
log.risk <- log(rep(c(rep(1, 70), rep(2, 30)),t))
x <- rnorm(n.all)
risk <- exp(log.risk + 0.1 * x)
mean <- E * risk
Y <- rpois(n=n.all, lambda=mean)
formula <- Y~ offset(log(E)) + x
     

#### Run the model     
model1 <- ST.knorrheld.main(formula, data=NULL, W=W, burnin=5000, n.sample=10000)

Run the code above in your browser using DataLab