gistfile1.txt

library(raster)
library(RandomFields)
library(R2OpenBUGS)
library(coda)
library(R2jags)
library(mvtnorm)

# Compare JAGS manual method to spatial.exp in BUGS

# simulate data
Res<-20

# Generate toy elevation data
Elevation<-raster(nrows=Res, ncols=Res, xmn=0, xmx=1, ymn=0, ymx=1, 
                  vals=NULL)
Elevation[]<-coordinates(Elevation)[,1]*coordinates(Elevation)[,2]

# Randomly select some pixels to act as observation coordinates
set.seed(1981)
Coords<-coordinates(Elevation)[sample(1:Res^2,100,replace=F),]

# Simulate a random field as error with 0 mean
lambda<-0.2
DistMat<-as.matrix(dist(coordinates(Elevation),diag=TRUE, upper=TRUE))
SIGMA <- exp(-DistMat/lambda)
mu <- rep(0, Res^2)

# sampling from the multivariate normal distribution
M <- Elevation
set.seed(1981)
M[] <- rmvnorm(1, mu, SIGMA)

# Simulate data using covariates and spatially structured error
alpha<- -2
beta<- 1
LogOdds<- alpha + beta*Elevation + M
Prob<-exp(LogOdds)/(1+exp(LogOdds))
DataProb<-extract(Prob,Coords)
Data<-data.frame(Coords=Coords,NPos=rbinom(100,100,DataProb))
Data$Elevation<-extract(Elevation,Coords)


# Prepare for JAGS and BUGs
SpatialTest_Data<- list(NClusters = nrow(Data),
                        Elevation = Data$Elevation,
                        y = Data$NPos,
                        Pop = rep(100,100),
                        Lat = Coords[,2],
                        Long = Coords[,1],
                        D = as.matrix(dist(Coords,diag=TRUE, upper=TRUE))) 
                        
                        
# Run JAGS model
cat("
    model
{
    # priors
    lambda ~ dunif(0.01, 100)
    alpha ~ dnorm(0, 0.001)
    beta ~ dnorm(0, 0.001)
    tau ~ dgamma(0.001,0.001)
    V ~ dgamma(0.001,0.001)
    
    for(i in 1:NClusters)
{
    # Vector of RE means
    mu[i]<-0
    
    #Define output as coming from binomial distribution
    y[i] ~ dbin(p[i], Pop[i])
    
    # Logistic model 
    logit(p[i]) <- alpha + beta[1]*Elevation[i] + u[i] 
}
    
    # derived quantities
    for(i in 1:NClusters)
{
    
    for(j in 1:NClusters)
{
    
    # turning the distance matrix to covariance matrix
    D.covar[i,j] <- V * exp(-D[i,j]/lambda)
}
}
    # turning covariances into precisions 
    D.tau[1:NClusters,1:NClusters] <- inverse(D.covar[1:NClusters,1:NClusters])
    
    # likelihood
    u[1:NClusters] ~ dmnorm(mu[], D.tau[,])
    }
    ", file="TEST_Spatial_Mod_JAGS.txt")


results_JAGS <-  jags(data=SpatialTest_Data, 
             parameters.to.save=c("lambda","V", "alpha","beta"),
             model.file="TEST_Spatial_Mod_JAGS.txt",
             n.iter=10000,
             n.chains=2,
             n.burnin=8000,
             n.thin=10,
             DIC=TRUE)

plot(as.mcmc(results_JAGS))
summary(as.mcmc(results_JAGS))