ortho.CAR<-function(O,E,x,s,np1,np2,
           spatial=T,ortho=F,
           sd.beta=1000,a.e=0.01,b.e=0.01,
           a.s=0.5,b.s=0.005,a.r=1,b.r=1,
           runs=5000,burn=1000){

    #DATA:
    #O = n x 1 matrix of counts
    #E = n x 1 matrix of expected counts (offsets)
    #x = n x p matrix of covariates
    #s = n x 1 matrix of spatial locations, s[i] \in {1,...,ns} 
    #The neighbor pairs are (np1[i],np2[i]),i=1,...,number of pairs

    #MODEL:
    #O ~ Pois(E*exp(x%*%beta + P%*%theta[s] + H))
    #P=diag(n)-x%*%solve(t(x)%*%x)%*%t(x)
    #H ~ N(0,taue)
    #theta~CAR(taus,rho)
    #beta~N(0,sd.beta)
    #taue~gamma(a.e,b.e)
    #taus~gamma(a.s,b.s)
    #rho~beta(a.r,b.r)

    #MODEL OPTIONS/HYPERPARAMETERS
    #spatial = F sets theta to zero
    #ortho = F sets P to diag(n)
    
    #Set up data:
    n<-nrow(x)
    p<-ncol(x)
    ns<-max(c(np1,np2,s))

    #spatial adjacency matrix:
    if(spatial){
     ADJs<-matrix(0,ns,ns)
     for(j in 1:length(np1)){
       ADJs[np1[j],np2[j]]<-ADJs[np2[j],np1[j]]<-1
     }
     Ms<-diag(apply(ADJs,2,sum))

     P<-matrix(0,n,ns)
     P[cbind(1:n,s)]<-1
     if(ortho){P<-(diag(n)-x%*%solve(t(x)%*%x)%*%t(x))%*%P}
    }

    #initial values   
    H<-log(O+1)-log(E+1)
    beta<-as.vector(lm(H~x-1)$coef)
    theta<-rep(0,ns)
    taue<-taus<-1
    rho<-0.95
   
    #keep track of stuff:
    keep.beta<-matrix(0,runs,p)
    theta.mn<-theta.var<-0*theta
    H.mn<-0*H
    keep.taue<-keep.taus<-keep.rho<-dev<-rep(0,runs)

    txx<-t(x)%*%x
    if(spatial){
     tPP<-t(P)%*%P
     M2<-diag(1/sqrt(diag(Ms)))
     ds<-eigen(M2%*%ADJs%*%M2)$values
     canrho<-detpart<-qbeta(seq(0.001,0.999,length=100),a.r,b.r)
     for(j in 1:length(canrho)){
      detpart[j]<-0.5*sum(log(1-canrho[j]*ds))
     }
     rm(M2,ds)
    }

    MH<-1
    #Start MCMC:
    for(i in 1:runs){

      ####################################################
      ######        Heterogeneity effects            #####
      ####################################################:
      mn<-x%*%beta
      if(spatial){mn<-mn+P%*%theta}
      sig<-1/sqrt(taue)
      canH<-rnorm(n,H,MH)
      R<-dpois(O,E*exp(canH),log=T)-dpois(O,E*exp(H),log=T)+
         dnorm(canH,mn,sig,log=T)-dnorm(H,mn,sig,log=T)
      acc<-runif(n)<exp(R)
      H<-ifelse(acc,canH,H)
      if(mean(acc)<0.35){MH<-MH*0.95}
      if(mean(acc)>0.45){MH<-MH*1.05}

      ####################################################
      ######         update regression coefs         #####
      ####################################################:
      V<-solve(taue*txx + diag(p)/sd.beta^2)
      if(spatial){M<-taue*t(x)%*%(H-P%*%theta)}
      if(!spatial){M<-taue*t(x)%*%H}
      beta<-V%*%M+t(chol(V))%*%rnorm(p)
      r<-H-x%*%beta

      ####################################################
      ######     update spatial random effects       #####
      ####################################################:

      if(spatial){
        V<-solve(taue*tPP + taus*(Ms-rho*ADJs))
        M<-taue*t(P)%*%r
        theta<-V%*%M+t(chol(V))%*%rnorm(ns)
        r<-r-P%*%theta
      }
      ####################################################
      ######        update the variances             #####
      ####################################################:
      SSE<-sum(r^2)
      taue<-rgamma(1,n/2+a.e,SSE/2+b.e)

      if(spatial){
       SS1<-t(theta)%*%Ms%*%theta
       SS2<-t(theta)%*%ADJs%*%theta
       taus<-rgamma(1,ns/2+a.s,(SS1-rho*SS2)/2+b.s)
       R<-detpart+0.5*taus*canrho*sum(SS2)
       rho<-sample(canrho,1,prob=exp(R-max(R)))
      }

      #keep track of stuff:
       keep.beta[i,]<-beta
       keep.taue[i]<-taue
       keep.taus[i]<-taus
       keep.rho[i]<-rho
       dev[i]<--2*sum(dpois(O,E*exp(H),log=T))

       if(i>burn){
          H.mn<-H.mn+H/(runs-burn)
          theta.mn<-theta.mn+theta/(runs-burn)
          theta.var<-theta.var+theta*theta/(runs-burn)
       }

      #Display current iteration:
      if(i%%1000==0){
        print(paste("Done with",i,"iterations"))
      }

   }

   theta.var<-theta.var-theta.mn^2
   dhat<--2*sum(dpois(O,E*exp(H.mn),log=T))
   dbar<-mean(dev[burn:runs])
   pD<-dbar-dhat
   DIC<-dbar+pD 

list(theta.mn=theta.mn,theta.var=theta.var,beta=keep.beta,
     taue=keep.taue,taus=keep.taus,rho=keep.rho,    
     dev=dev,dbar=dbar,pD=pD,DIC=DIC)}    

#OUTPUT:
#theta.mn is the posterior mean of the CAR random effects
#theta.var is the posterior variance of the CAR random effects
#taue is the draws of taue
#taus is the draws of taus
#rho is the draws of rho
#beta is the draws of beta
#dev is the draws of the deviance
#DIC, dbar, and pD are DIC statistics