source("C:\\...\\orthoCARpoisson.R")

n<-100        #number of observations:
ns<-100       #number of counties
p<-3          #number of covariates
runs<-25000   #number of MCMC iterations
burn<-5000    #length of burn-in


#True values
theta<-sin(2*pi*5*1:ns/ns)
beta<-1:p/p/10
E<-rep(10,n)



#Generate data
s<-1:ns
x<-matrix(rnorm(n*p),n,p)
x[,1]<-1  #intercept:
x[,p]<-rnorm(n,theta[s],0.5) #correlated with the spatial trend
O<-rpois(n,E*exp(x%*%beta+theta[s]))



#spatial grid is a line
np1<-2:ns-1
np2<-2:ns


#fit three models
fit1<-ortho.CAR(O,E,x,s,np1,np2,runs=runs,burn=burn,spatial=F)
fit2<-ortho.CAR(O,E,x,s,np1,np2,runs=runs,burn=burn,spatial=T,ortho=F)
fit3<-ortho.CAR(O,E,x,s,np1,np2,runs=runs,burn=burn,spatial=T,ortho=T)

print(summary(glm(O~x-1,family="poisson",offset=log(E))))


#Compare the models with DIC:
print(" ");print(" ");print(" ")
print("***** NON-SPATIAL MODEL ****")
print("D-bar, pD, DIC")
print(round(c(fit1$dbar,fit1$pD,fit1$DIC)))
print("95% intervals for beta")
print(round(apply(fit1$beta[burn:runs,],2,quantile,c(0.025,0.5,0.975)),3))


print(" ");print(" ");print(" ")
print("***** SPATIAL MODEL ****")
print("D-bar, pD, DIC")
print(round(c(fit2$dbar,fit2$pD,fit2$DIC)))
print("95% intervals for beta")
print(round(apply(fit2$beta[burn:runs,],2,quantile,c(0.025,0.5,0.975)),3))


print(" ");print(" ");print(" ")
print("***** ORTHOGONALIZED SPATIAL MODEL ****")
print("D-bar, pD, DIC")
print(round(c(fit3$dbar,fit3$pD,fit3$DIC)))
print("95% intervals for beta")
print(round(apply(fit3$beta[burn:runs,],2,quantile,c(0.025,0.5,0.975)),3))