model{
#Tree survival model
 for(i in 1:trees){
   tree.term[i]~dnorm(0, tau.tree)    
   pr.alive[i,1]<-1
   for(t in 2:last[i]){
                 rabdens[i,t]<-(warrendat[site.id[i], t-1]+warrendat[site.id[i], t])/2 
                 logit(pr.surv[i,t])<-
                 beta[2]*(height[i,t-1]/1000)+            #effect of previous height on survival
                 beta[3]*DS[i]+                        	#effect of being a drooping sheoak
                 beta[4]*GW[i]+                         #effect of being a golden wattle
                 beta[5]*unguard[i]+                      #effect of being guarded
  beta[6]*rabdens[i,t]+    		#effect of rabbit density at t-1                       
  beta[7]*unguard[i]*rabdens[i,t]+            #guard by ripping interaction   
                 beta[8]*(height[i,t-1]/1000)*rabdens[i,t]+ #height by density interaction
                 beta[9]*(height[i,t-1]/1000)*unguard[i]+  #height x guard interaction
                 beta[10]*(height[i,t-1]/1000)*rabdens[i,t]*guard[i] + #height x guard x density interaction 
                 site.term[site.id[i]]  +                              #site level random effects  
                 cluster.term[cluster.id[i]]      +                   #cluster level random effects 
                 tree.term[i]    + #individual-level random effects
                 time.term[t-1]                   #time-step level random effect                 
                 pr.surv.corrected[i,t]<-pow(pr.surv[i,t], t.ints[t-1])
                 Y[i,t]~dbern(pr.surv.corrected[i,t])

#POSTERIOR PREDICTIVE CHECK - compute chi-sq discrepancy for each tree
                chi.real[i,t]<-pow(Y[i,t]-pr.surv.corrected[i,t], 2)/pr.surv.corrected[i,t]
                 Y.fake[i,t]~dbern(pr.surv.corrected[i,t])
                 chi.fake[i,t]<-pow(Y.fake[i,t]-pr.surv.corrected[i,t], 2)/pr.surv.corrected[i,t]
                 
  }
#Above-ground biomass model (for surviving trees)
logBiomass[i]~dnorm(mu[i], tau.biomass)
         mu[i]<-raneff[site.id[i]]+
         BB[2] *     DS[i]+
         BB[3] *     GW[i]+
         BB[4] *  unguard[i]+
         BB[5] * mean(warrendat[site.id[i],])+         
         BB[6] * mean(warrendat[site.id[i],])*unguard[i]

#POSTERIOR PREDICTIVE CHECK - compute squared residuals for real and fake biomass data
   logBiomass.fake[i]~dnorm(mu[i], tau.biomass)
   sq.residual.real[i]<-pow(mu[i]-logBiomass[i], 2)
   sq.residual.fake[i]<-pow(mu[i]-logBiomass.fake[i], 2)
  }

#PRIORS for tree survival model
#cluster level variation in survival probability (#hierarchical centering applied)
for(c in 1:clusters){ cluster.term[c]~dnorm(0, tau.clus) }  
  tau.clus<-1/pow(sigma.clus, 2)
  sigma.clus~dunif(0, 100)

#site level variation in survival probability  (#hierarchical centering not applied)
for(s in 1:sites){ site.term[s]~dnorm(beta[1], tau.site) } 
  tau.site<-1/pow(sigma.site, 2)
  sigma.site~dunif(0, 100)

for(aa in 1:trees){ tree.term[aa]~dnorm(0, tau.tree) } 
  tau.tree<-1/pow(sigma.tree, 2)
  sigma.tree~dunif(0, 100)

for(tt in 1:7) {time.term[tt]~dnorm(0, tau.time)}
   tau.time<-1/pow(sigma.time, 2)
   sigma.time~dunif(0, 100)  

#regression parameters
for(b in 1:10){beta[b]~dnorm(0, 0.0001) }

#PRIORS for above-ground biomass model (for surviving trees)
#regression parameters
for(k in 1:6) {BB[k]~dnorm(0, 0.0001)}

for(ss in 1:sites){raneff[ss]~dnorm(BB[1], tau.site.biomass) } #site level random effects on biomass

#site/biomass error term
tau.site.biomass<-1/pow(sigma.site.biomass, 2)
sigma.site.biomass~dunif(0, 100)

#residual biomass error term
tau.biomass<-1/pow(sigma.biomass, 2)
sigma.biomass~dunif(0, 100)

#compute the Bayesian p-value to check goodness of fit for tree survival model
for(kkk in 1:trees){
rowsums_chi.fake[kkk]<-sum(chi.fake[kkk, 2:last[kkk]])
rowsums_chi.real[kkk]<-sum(chi.real[kkk, 2:last[kkk]])
}

chisq.fake<-sum(rowsums_chi.fake)
chisq.real<-sum(rowsums_chi.real)
#will be one if chisq. discrepancy of fake data more than real
p.val<-step(chisq.fake-chisq.real)

#compute the Bayesian p-value to check goodness of fit for above-ground biomass model
SS.biomass.real<-sum(sq.residual.real)
SS.biomass.fake<-sum(sq.residual.fake)
#will be one if resid. variance of fake data more than real
p.val.biomass<-step(SS.biomass.fake-SS.biomass.real) 

}