ortho.CAR<-function(Y,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:
    #Y = n x 1 matrix of observations
    #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:
    #Y ~ N(x%*%beta + P%*%theta[s], taue*diag(n))
    #P=diag(n)-x%*%solve(t(x)%*%x)%*%t(x)
    #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   
    beta<-as.vector(lm(Y~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
    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)
    }

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

      ####################################################
      ######         update regression coefs         #####
      ####################################################:
      V<-solve(taue*txx + diag(p)/sd.beta^2)
      if(spatial){M<-taue*t(x)%*%(Y-P%*%theta)}
      if(!spatial){M<-taue*t(x)%*%Y}
      beta<-V%*%M+t(chol(V))%*%rnorm(p)
      r<-Y-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(dnorm(r,0,1/sqrt(taue),log=T))

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

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

   }

   theta.var<-theta.var-theta.mn^2
   beta.mn<-apply(keep.beta[burn:runs,],2,mean)
   sigma<-mean(1/sqrt(keep.taue[burn:runs]))
   r<-Y-x%*%beta.mn
   if(spatial){r<-r-P%*%theta.mn}
   dhat<-sum(-2*dnorm(r,0,sigma,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