#
#  R Script to fit Occupancy model to Willow Tit Data
#  
#  Model Selection Approach: Information Criteria 
#


###
###  Load Covariates and Occupancy Data
###

wt.df=read.csv("wt.csv",header=TRUE)
head(wt.df)

n=200
Y=as.matrix(wt.df[1:n,1:3])
y=apply(wt.df[1:n,1:3],1,sum,na.rm=TRUE)
J=apply(!is.na(wt.df[1:n,1:3]),1,sum)

X=matrix(1,n,2)
X[,1]=scale(wt.df$elev[1:n])
X[,2]=scale(wt.df$forest[1:n])

###
###  Fit Each Probit Occupancy Model in C-V Loop 
###

n.mcmc=160000
s2.beta=1.5^2

WAIC.vec=rep(0,4)
DIC.vec=rep(0,4)
D.vec=rep(0,4)
source("occ.aux.mcmc.R")

tmp.time=proc.time()
tmp.out.1=occ.aux.mcmc(Y,J,X,1,1,1,n.mcmc,s2.beta,null.model=TRUE,no.print=TRUE,get.D=TRUE,no.pred=TRUE)   
WAIC.vec[1]=tmp.out.1$WAIC.2
DIC.vec[1]=tmp.out.1$DIC
D.vec[1]=tmp.out.1$D
tmp.out.2=occ.aux.mcmc(Y,J,X[,1],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,get.D=TRUE,no.pred=TRUE)
WAIC.vec[2]=tmp.out.2$WAIC.2
DIC.vec[2]=tmp.out.2$DIC
D.vec[2]=tmp.out.2$D
tmp.out.3=occ.aux.mcmc(Y,J,X[,2],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,get.D=TRUE,no.pred=TRUE)   
WAIC.vec[3]=tmp.out.3$WAIC.2
DIC.vec[3]=tmp.out.3$DIC
D.vec[3]=tmp.out.3$D
tmp.out.4=occ.aux.mcmc(Y,J,X[,1:2],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,get.D=TRUE,no.pred=TRUE)   
WAIC.vec[4]=tmp.out.4$WAIC.2
DIC.vec[4]=tmp.out.4$DIC
D.vec[4]=tmp.out.4$D
time.1=proc.time()-tmp.time
time.1

WAIC.vec
DIC.vec
D.vec