###Part I: Preparing the data
rm(list=ls())
  library(nlme)
  library(bbmle)
  library(lme4)     # for fitting GLMMs
  library(lattice)  # for the xyplot function

##Set your working directory manually or using the following commands## 
#getwd()
#setwd("pathway")
#setwd("C:/Users/Pedro/Documents/Classes/Methods in Ecology II/Demos/")

## read the Hypericum data from file

## read the Hypericum data from file
orig_data <- read.table("hypericum_data_94_07.txt", header=T)
dt <- subset(orig_data, !is.na(ht_init) & !is.na(st_init) & rp_init > 0 & year<1997 )
yr <- unique(dt$year)
dt$lgh <- log(dt$ht_init)
dt$lfr <- log(dt$rp_init)
dt$stems <- dt$st_init
site <- unique(dt$bald)
table(dt$bald,dt$fire_year)
dt$TSF <- 1
dt$TSF[dt$fire_year <1987] <-2 
dt$TSF[dt$fire_year <1973] <-3 
dt$TSF <- factor(dt$TSF)
dt$fyear <- factor(dt$year)
dt$fbald <- factor(dt$bald)
dt$surv <-1
dt$surv[dt$fate =="rip"] <- 0
table(dt$surv,dt$fate)
I <- order(dt$lgh)
lgh <- sort(dt$lgh)

table(dt$bald,dt$TSF)
tsf <- unique(dt$TSF)
tsf <- sort(tsf)
TSF <-dt$TS
boxplot(dt$lgh~dt$TSF)
summary(lm(dt$lgh~factor(dt$TSF)))
pairs(subset(dt,select=c(ht_init, rp_init,stems)))
####################################################################
dt$stems[dt$stems>8] <- 8
dt$stems <- factor(dt$stems)
dt$lghc <- dt$lgh- mean(dt$lgh)

######################  RANDOM MODELS ##############################################
require(optimx)
m1 <- glm(surv~lghc*TSF*stems,data=dt,family =binomial)
m2 <- glmer(surv~lghc*TSF*stems + (1|fbald) + (1|year),data=dt,family =binomial)
m2_nlminb <- update(m2,control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
m3 <- glmer(surv~lghc*TSF*stems + (lghc|fbald)+ (lghc|year),data=dt,family =binomial)
m3_nlminb <- update(m3,control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))

AICtab(m1,m2_nlminb,m3_nlminb,weights=TRUE,base = TRUE)


############################# FIXED MODELS ###########################################

M11 <- glmer(surv~lghc*TSF*stems + (1|fbald) + (1|year),data=dt,family=binomial)
M11_nlminb <- update(M11,control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
M13 <- glmer(surv~lghc+TSF*stems + (1|fbald) + (1|year),data=dt,family =binomial)
M13_nlminb <- update(M13,control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
M14 <- glmer(surv~lghc+TSF+stems + (1|fbald) + (1|year),data=dt,family =binomial)
M15 <- glmer(surv~lghc+TSF + (1|fbald) + (1|year),data=dt,family =binomial)
M16 <- glmer(surv~lghc+stems + (1|fbald) + (1|year),data=dt,family =binomial)
M17 <- glmer(surv~lghc*stems + (1|fbald) + (1|year),data=dt,family =binomial)
M18 <- glmer(surv~lghc*stems + TSF + (1|fbald) + (1|year),data=dt,family =binomial)
M19 <- glmer(surv~lghc*TSF + stems + (1|fbald) + (1|year),data=dt,family =binomial)
M20 <- glmer(surv~lghc*TSF + stems*TSF + (1|fbald) + (1|year),data=dt,family =binomial)
M20_nlminb <- update(M20,control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))

AICtab(M11_nlminb,M13_nlminb,M14,M15,M16,M17,M18,M19,M20_nlminb,weights=TRUE,base = TRUE)

summary(M20_nlminb)


#######################################################################
color_eren <- c("yellow2","yellowgreen","forestgreen","turquoise3","dodgerblue3",
                "darkslateblue","blueviolet","red4")
par(mfrow=c (1,1)) 
plot(rep(1,8),col=color_eren,pch=19,cex=3)
min(dt$lghc)
max(dt$lghc)
x<- seq(-1.4, 0.85, 0.01)
par(mfrow=c (1,3),bg="lightgray")

for (i in 1:3){
plot(exp(dt$lghc+mean(dt$lgh)),dt$surv,type="n",ylim=c(0,1),xlim=exp(c(-1.4, 0.85)+mean(dt$lgh)),main=paste("TSF",tsf[i]),ylab="P(survival)",xlab="height(cm)")
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = 
"grey")
for (j in 1:8){
 coeff <- fixef(M20_nlminb)
 if (i == 1) {
      coeff_tsf <- 0
      coeff_TSFxstems <- 0
      coeff_tsfxlgh <- 0
      }
 if (i == 2) {
      coeff_tsf <- coeff[3]
      if (j == 1) {coeff_TSFxstems <- 0}
      coeff_tsfxlgh <- coeff[12]
      if (j == 1) {coeff_TSFxstems <- 0}
      if (j == 2) {coeff_TSFxstems <- coeff[14]}
      if (j == 3) {coeff_TSFxstems <- coeff[16]}
 	if (j == 4) {coeff_TSFxstems <- coeff[18]}
      if (j == 5) {coeff_TSFxstems <- coeff[20]}
      if (j == 6) {coeff_TSFxstems <- coeff[22]}
 	if (j == 7) {coeff_TSFxstems <- coeff[24]}
      if (j == 8) {coeff_TSFxstems <- coeff[26]}
      }
 if (i == 3) {
      coeff_tsf <- coeff[4]
      if (j == 1) {coeff_TSFxstems <- 0}
      coeff_tsfxlgh <- coeff[13]
 	if (j == 2) {coeff_TSFxstems <- coeff[15]}
      if (j == 3) {coeff_TSFxstems <- coeff[17]}
 	if (j == 4) {coeff_TSFxstems <- coeff[19]}
      if (j == 5) {coeff_TSFxstems <- coeff[21]}
      if (j == 6) {coeff_TSFxstems <- coeff[23]}
 	if (j == 7) {coeff_TSFxstems <- coeff[25]}
      if (j == 8) {coeff_TSFxstems <- coeff[27]}
      }
 if (j == 1) {coeff_stems <- 0}
 if (j == 2) {coeff_stems <- coeff[5]}
 if (j == 3) {coeff_stems <- coeff[6]}
 if (j == 4) {coeff_stems <- coeff[7]}
 if (j == 5) {coeff_stems <- coeff[8]}
 if (j == 6) {coeff_stems <- coeff[9]}
 if (j == 7) {coeff_stems <- coeff[10]}
 if (j == 8) {coeff_stems <- coeff[11]}
        odds_ratio <- coeff[1]+coeff[2]*x+coeff_tsf+coeff_stems+coeff_tsfxlgh*x+coeff_TSFxstems
      prob <- 1/(1 + (1/exp(odds_ratio)))
lines(exp(x+mean(dt$lgh)), prob,col=color_eren[j],lwd=3)
}}

###########################################################
par(mfrow=c (1,4))
hist(residuals(M20_nlminb),main ="residuals")
plot(dt$lgh,residuals(M20_nlminb))
plot(dt$stems,residuals(M20_nlminb))
plot(TSF,residuals(M20_nlminb))


#######################################################################


##### Bayesian Mixed Model #############
###  difuse priors for fixed effects and for random effects ##
library(jagsUI) 

dt$ppn[dt$bald==1] <- 1
dt$ppn[dt$bald==29] <- 2
dt$ppn[dt$bald==32] <- 3
dt$ppn[dt$bald==42] <- 4
dt$ppn[dt$bald==50] <- 5
dt$ppn[dt$bald==57] <- 6
dt$ppn[dt$bald==59] <- 7
dt$ppn[dt$bald==62] <- 8
dt$ppn[dt$bald==67] <- 9
dt$ppn[dt$bald==87] <- 10
dt$ppn[dt$bald==88] <- 11
dt$ppn[dt$bald==91] <- 12
dt$ppn[dt$bald==93] <- 13
dt$ppn[dt$bald==103] <-14

dt$yrn[dt$year==1994] <- 1
dt$yrn[dt$year==1995] <- 2
dt$yrn[dt$year==1996] <- 3

X <- model.matrix(~surv,data=dt)
K <- ncol(X)
Npop <- length(unique(dt$bald))
Nyear <- length(unique(dt$year))
Nstems <- length(unique(dt$stems))
Ntsf <- length(unique(dt$TSF))
n <- nrow(dt)


#################################################################

################# intercept

model =  "Model_w_year binary intercept.R"

# Bundle data
win.data1=list(
  Y = dt$surv,
  X = dt$lghc,
  N = n,
  pp = as.numeric(dt$ppn),
  yy = as.numeric(dt$yrn),
  tsf = as.numeric(dt$TSF),
  stem = as.numeric(dt$stems),
  Npop = Npop,
  Nyear = Nyear,
  Nstems= Nstems,
  Ntsf = Ntsf
)

# Inits function 
inits <- function(){
  list(sigma.plot.a = runif(1,0.001,10),
       sigma.year.a = runif(1,0.001,10),
       sigma.eps = runif(1,0.001,10)      
  )}


# Parameters to estimate
params <- c("Betas1","Betas2","Betas3","Betas4","Betas5","Betas6","alpha.p","alpha.y",
            "sigma.plot.a","sigma.year.a","sigma.eps")

#jm=jags(model,data=win.data,inits,n.chains=3,n.adapt=5000) 
jm=jags(win.data1,inits,params,model.file=model,n.chains=3, n.iter=10000, n.burnin=5000, n.thin=100)

jm <- update(jm, n.iter=100000)
jm								                      #Call model results
#plot(jm)
#################################################################