##################################################
#### Run the model on simulated data on a lattice
##################################################
#### Set up a square lattice region
x.easting <- 1:10
x.northing <- 1:10
Grid <- expand.grid(x.easting, x.northing)
n <- nrow(Grid)
#### 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 the covariates and response data
x1 <- rnorm(n)
x2 <- rnorm(n)
theta <- rnorm(n, sd=0.05)
phi <- mvrnorm(n=1, mu=rep(0,n), Sigma=0.4 * exp(-0.1 * distance))
fitted <- -0.2 + 0.1 * x1 + 0.1*x2 + theta + phi
Y <- rnorm(n=n, mean=fitted, sd=rep(1,n))
#### Run the BYM model
#### Let the function randomly generate starting values for the parameters
#### Use the default priors specified by the function (for details see the help files)
formula <- Y ~ x1 + x2
model <- gaussian.bymCAR(formula=formula, W=W, burnin=5000, n.sample=10000)
model <- gaussian.bymCAR(formula=formula, W=W, burnin=20, n.sample=50)
Run the code above in your browser using DataLab