library(plyr)
require(gdata)
library(R2WinBUGS)

##record of deaths
##read in data
deaths<-read.csv("\\\\psf/Home/Documents/Bird_Facade_Data/deaths.csv",head=T)
deaths<-deaths[,2:6]
colnames(deaths)<-c("Year","Month","Building","Facade","Species") ##don't care about species
##information on number of buildings and number of facades (## of observations)
windows<-read.csv("Y:/Documents/Bird_Facade_Data/windows.csv",head=T)
colnames(windows)<-c("Building","Facade","SqM")

##fill in NA
for(i in 1:dim(windows)[1]){ ##each row
	windows[i,1]<-ifelse(is.na(windows[i,1])==T,windows[i-1,1],windows[i,1])
}

##begins October 2013 - May 2015 -- total of 3+12+5 months==20
##20 time points

deaths$UniqueTime<-paste(deaths$Month,deaths$Year,sep="")
deaths$UniqueTimeBuilding<-paste(deaths$UniqueTime,"B",deaths$Building,sep="")
deaths$UniqueTimeFacade<-paste(deaths$UniqueTime,"B",deaths$Building,"F",deaths$Facade,sep="")

##want table of times (20), for each building, for each --- 100 facades * 20 time points = 2000 points
##let's do N x T where N=100 facades and T=20 points

####building dat
datb<-matrix("NA",ncol=length(unique(deaths$UniqueTime)),nrow=length(unique(windows$Building)))
##counts number of bird deaths per building for each time point (month in year)
for(j in 1:length(unique(windows$Building))){ ##N=22
	for(i in 1:length(unique(deaths$UniqueTime))){ ##T=20 		
	datb[j,i]<-ifelse((paste(unique(deaths$UniqueTime)[i],"B",unique(windows$Building)[j],sep="")%in%deaths$UniqueTimeBuilding)==TRUE,
	
					count(deaths,("UniqueTimeBuilding"))[which(paste(unique(deaths$UniqueTime)[i],"B",unique(windows$Building)[j],sep="")==count(deaths,("UniqueTimeBuilding"))[,1]),2],0)
	}
}
datb<-as.data.frame(datb)
colnames(datb)<-unique(deaths$UniqueTime)

##read in covariates by building
covariates<-read.csv("Y:/Documents/Bird_Facade_Data/VTCRC_collisions.csv",header=T)

datb$BuildingNumber<-as.character(unique(windows$Building))
datb<-merge(datb,covariates,by="BuildingNumber")


###figure out number of observations per building per time point
##sorted sheet 5 by date
visits<-read.csv("Y:/Documents/Bird_Facade_Data/Schneider_et al_data_R5_for_Keith.csv",head=T)
visits<-visits[,c(2:5,9)]
visits$Month<-toupper(visits$Month)
visits$UniqueTime<-paste(substr(visits$Month,1,3),visits$Year,sep="")
visits$UniqueTime[which(visits$UniqueTime=="JUL204")]<-"JUL2014"
visits<-visits[-which(visits$UniqueTime=="NA"),] ##get rid of observations with no time
visits<-visits[-which(visits$Building=="x"),] ##get rid of observations where we don't know the building
##setdiff(as.vector(unique(visits$Building)),as.vector(unique(deaths$Building)))
##buildings 5, 19, and 25 aren't in deaths set
visits<-visits[-which(visits$Building%in%c(5,19,25)),] ##get rid of buildings 5, 19, and 25 since we have no info on them
visits$Building<-factor(visits$Building)
visits$UniqueTime<-as.factor(visits$UniqueTime)
levels(visits$UniqueTime)<-c("OCT2013",           "NOV2013",           "DEC2013",           "JAN2014",          
"FEB2014",           "MAR2014",           "APR2014",           "MAY2014",           "JUN2014",          
"JUL2014",           "AUG2014",           "SEPT2014",          "OCT2014",           "NOV2014",          
"DEC2014",           "JAN2015",           "FEB2015",           "MAR2015",           "APR2015",          
"MAY2015")

##time is i; building is j
obs<-matrix("NA",ncol=length(unique(visits$UniqueTime)),nrow=length(unique(visits$Building)))
##counts number of bird deaths per building for each time point (month in year)
for(j in 1:length(levels(visits$Building))){ ##N=22
	for(i in 1:length(levels(visits$UniqueTime))){ ##T=20 		
	obs[j,i]<-ifelse(levels(visits$UniqueTime)[i]%in%count(visits,c("UniqueTime","Building"))[,1]&levels(visits$Building)[j]%in%count(visits,c("UniqueTime","Building"))[,2],###if it's in the table
		count(visits,c("UniqueTime","Building"))[which(count(visits,c("UniqueTime","Building"))[,1]==levels(visits$UniqueTime)[i]&count(visits,c("UniqueTime","Building"))[,2]==levels(visits$Building)[j]),3],"0")
	}
}

obs<-as.data.frame(obs)
##change factor columns to numeric

for(i in 1:ncol(obs)) {
    obs[,i] <- as.numeric(as.character(obs[,i]))
}

obs[is.na(obs)]<-0 ##no surveys done for these months at these buildings
colnames(obs)<-levels(visits$UniqueTime)
obs$BuildingNumber<-levels(visits$Building)

##for BRugs
T<-length(unique(deaths$UniqueTime)) ##number of time points
N<-dim(datb)[1] ##number of buildings
observed<-datb[,2:21] ##drop ID column
for(i in 1:ncol(observed)) {
    observed[,i] <- as.numeric(as.character(observed[,i]))
}
observed<-as.matrix(observed)
pop<-obs[,1:T] ##population as in total number of surveys for each month
sum.count<-apply(observed,2,sum) ##sums total number of birds death for each month
sum.count.pop<-apply(pop,2,sum) ##sums total number of surveys for each month
r<-sum.count/sum.count.pop ##rate of collision/death per survey for each month
expected<-as.matrix(pop)*r ##N*T %*% T
###add more covariates

##ZIP model
d<-list(N=N, T=T,
	observed=observed,
	expected=expected,
	x1=datb$PercentTree,
	x2=datb$PercentGrass,
	x3=as.vector(scale(datb$WindowArea))	
)

baymodel<-function(){
for(i in 1:N) {
	for( t in 1:T){
		observed[i,t] ~ dpois(mu[i,t])
		log(theta[i,t]) <- log(expected[i,t])+beta1*x1[i]+beta2*x2[i]+beta3*x3[i]+delta[t]
		mu[i,t] <- (1-W[i,t])*theta[i,t]+0.00001
		W[i,t]~dbern(p[i])
	}
	logit(p[i])<-gamma[i]
	gamma[i]<-beta1g*x1[i]+beta2g*x2[i]+beta3g*x3[i]
	}
##poisson parameters
beta1~dnorm(0,1.0E-2)
beta2~dnorm(0,1.0E-2)
beta3~dnorm(0,1.0E-2)

#binomial paramters
beta1g~dnorm(0,1.0E-2)
beta2g~dnorm(0,1.0E-2)
beta3g~dnorm(0,1.0E-2)

mu.delta~dnorm(0,1)
tau.delta<-0.001#dgamma(0.02,0.02)
pred.delta[1]<-mu.delta
delta[1]~dnorm(pred.delta[1],tau.delta)

for(t in 2:T){
pred.delta[t]<-delta[t-1]
delta[t]~dnorm(pred.delta[t],tau.delta)
}
}

modelfile <- "Y:/Documents/Bird_Facade_Data/model.txt"
write.model(baymodel, modelfile)

init1 <- list(beta1=0,beta2=0,beta3=0,delta=rep(0,T),#p=rep(0.5,N))#,
	 beta1g=0,beta2g=0,beta3g=0)
inits<-list(init1)

wdir <- "Y:/Documents/Bird_Facade_Data/Dir"
if(!file.exists(wdir)){dir.create(wdir)}
BugsDir <- "\\\\psf/Home/Applications/winbugs14_unrestricted/WinBUGS14"

bugs1 <- bugs(data = d, inits = inits,
   working.directory = wdir,
   parameters.to.save = c("theta","beta1","beta2","beta3","delta","p","beta1g","beta2g","beta3g"),
   n.chains = 1, n.iter = 30000, n.burnin = 15000, n.thin = 5,
   model.file = modelfile,
   bugs.directory = BugsDir,debug=F)


bugs1$DIC
output1 <- read.coda("Y:/Documents/Bird_Facade_Data/Dir/coda1.txt",
		"Y:/Documents/Bird_Facade_Data/Dir/codaIndex.txt")
mcmc <- mcmc.list(output1)
plot(mcmc[,c("beta1","beta2","beta3")])
plot(mcmc[,c("beta1g","beta2g","beta3g")])
plot(mcmc[,c("delta[1]","delta[3]","delta[7]")])
plot(mcmc[,c("p[1]","p[10]","p[20]")])
geweke.diag(mcmc[,c("delta[1]","delta[3]","beta1","beta2","beta3","p[1]","p[10]","p[20]","beta1g","beta2g","beta3g")])

exp(bugs1$mean$beta1);quantile(exp(bugs1$sims.list$beta1),c(0.025,.975))
exp(bugs1$mean$beta2);quantile(exp(bugs1$sims.list$beta2),c(0.025,.975))
exp(bugs1$mean$beta3);quantile(exp(bugs1$sims.list$beta3),c(0.025,.975))

scale(datb$WindowArea)
##https://stats.stackexchange.com/questions/254934/what-is-the-interpretation-of-scaled-regression-coefficients-when-only-the-predi


exp(bugs1$mean$beta1g);quantile(exp(bugs1$sims.list$beta1g),c(0.025,.975))
exp(bugs1$mean$beta2g);quantile(exp(bugs1$sims.list$beta2g),c(0.025,.975))
exp(bugs1$mean$beta3g);quantile(exp(bugs1$sims.list$beta3g),c(0.025,.975))