# Model type: RDHuggins 
# Parameters: S(survival),GammaDoublePrime(emigration),GammaPrime(return),p(captureprobability), c(recapture probability)

library(RMark)
library(ggplot2)
library(knitr)
library(ggplot2)

# Stage 1 of the Sequential Modelling Process ------------------------------
# Importing Data
dater<-import.chdata("MODIFYDIRECTORY/inp.txt",header=TRUE, use.comments = TRUE)
# Area 1 removed due to small sample sizes causing model fit problems.
dater<-subset(dater,AREA!=1)
dater<-subset(dater,ch!="000000000000000000000000000000000000000000000000000000000000000000000000")
intervals<-c(rep(c(0,0,0,0,0,0,0,0,0,0,0,1),5),c(0,0,0,0,0,0,0,0,0,0,0))

data.processed=process.data(dater,model="RDHuggins",time.intervals=intervals,groups=c("AREA","GRID","Prebait"))
data.processed$group.covariates$AREA<-relevel(data.processed$group.covariates$AREA,ref='3')


# Parameter Fix Function --------------------------------------------------


fix.params=function(data.ddl)
{
  # Fix survival estimates for area 1 and 2 from the first three years
  data.ddl$S$fix=NA
  data.ddl$S$fix[data.ddl$S$AREA==1 & data.ddl$S$time%in%1:3]=0
  data.ddl$S$fix[data.ddl$S$AREA==2 & data.ddl$S$time%in%1:3]=0
  data.ddl$S$time<-relevel(data.ddl$S$time,ref='5')
  
  data.ddl$GammaDoublePrime$fix=NA
  data.ddl$GammaDoublePrime$fix[data.ddl$GammaDoublePrime$AREA==1 & data.ddl$GammaDoublePrime$time%in%1:4]=0
  data.ddl$GammaDoublePrime$fix[data.ddl$GammaDoublePrime$AREA==2 & data.ddl$GammaDoublePrime$time%in%1:4]=0
  
  data.ddl$GammaPrime$fix=NA
  data.ddl$GammaPrime$fix[data.ddl$GammaPrime$AREA==1 & data.ddl$GammaPrime$time%in%1:5]=0
  data.ddl$GammaPrime$fix[data.ddl$GammaPrime$AREA==2 & data.ddl$GammaPrime$time%in%1:5]=0
  
  data.ddl$p$fix=NA
  data.ddl$p$fix[data.ddl$p$AREA==1 & data.ddl$p$session%in%1:3]=0
  data.ddl$p$fix[data.ddl$p$AREA==2 & data.ddl$p$session%in%1:3]=0
  data.ddl$p$fix[data.ddl$p$AREA==3 & data.ddl$p$session==1 & data.ddl$p$time == c(1,2,3,4,5,6,7,8)]=0
  data.ddl$p$session<-relevel(data.ddl$p$session,ref='6')
  
  data.ddl$c$fix=NA
  data.ddl$c$fix[data.ddl$c$AREA==1 & data.ddl$c$session%in%1:3]=0
  data.ddl$c$fix[data.ddl$c$AREA==2 & data.ddl$c$session%in%1:3]=0
  data.ddl$c$session<-relevel(data.ddl$c$session,ref='6')
  return(data.ddl)
}

# Step 1: Global ----------------------------------------------------------


# Step 1 
# Fit global model 

global.model=function(data.processed)
{
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  #
  # Define parameter models
  #
  c.global=list(formula=~GRID*Time)
  p.global=list(formula=~GRID*Time)
  S.global=list(formula=~GRID*time)
  Gamma.Prime=list(formula=~1)
  Gamma.Double.Prime=list(formula=~1)
  # Run global model
  mod_global<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime,GammaPrime=Gamma.Prime))
  
  
  return(collect.models())
  
}

global.results=global.model(data.processed)

(global.table<-kable(model.table(model.list = global.results, sort = TRUE, adjust = TRUE,
                                 pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                 model.name = TRUE)))


# Step 2: E and I  --------------------------------------------------------


# Step 2
# Fit gamma models
Gamma.models=function(data.processed)
{
  
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  #
  # Define parameter models
  #
  c.global=list(formula=~GRID*Time)
  p.global=list(formula=~GRID*Time)
  S.global=list(formula=~GRID*time)
  
  Gamma.Prime=list(formula=~1)
  Gamma.Double.Prime=list(formula=~1)
  Gammas.Zero=list(formula=~1,fixed=0)
  Gammas.Double.Zero=list(formula=~1,fixed=0)
  Gammas.random=list(formula=~1, share=TRUE)
  Gamma.Prime.GRID=list(formula=~GRID)
  Gamma.Double.Prime.GRID=list(formula=~GRID)
  Gammas.shared.GRID=list(formula=~GRID, share=TRUE)
  
  # Candidate Models
  mod_separate.gammas<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime,GammaPrime=Gamma.Prime))
  mod_gammas.random<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gammas.random))
  mod_gammas.zero<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gammas.Double.Zero,GammaPrime=Gammas.Zero))
  mod_Gamma.Double.GRID<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.GRID,GammaPrime=Gamma.Prime))
  mod_Gamma.GRID<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime,GammaPrime=Gamma.Prime.GRID))
  mod_Gamma.Shared.GRID<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gammas.shared.GRID))
  mod_Gammas.GRID<-mark(data.processed,data.ddl,model.parameters=list(c=c.global, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.GRID,GammaPrime=Gamma.Prime.GRID))
  
  return(collect.models())
  
}
gamma.results=Gamma.models(data.processed)

(gamma.table<-kable(model.table(model.list = gamma.results, sort = TRUE, adjust = TRUE,
                                pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                model.name = FALSE)))

# Step 3: Behavioral ------------------------------------------------------


# Step 3
# Fit behavior models
behavior.models=function(data.processed)
{
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  #
  # Define parameter models
  #
  
  S.global=list(formula=~GRID*time)
  Gamma.Prime.best=list(formula=~1,fixed=0)
  Gamma.Double.Prime.best=list(formula=~1,fixed=0)
  
  c.dot=list(formula=~1)
  p.dot=list(formula=~1)
  c.Time=list(formula=~Time)
  p.Time=list(formula=~Time)
  p.c.share=list(formula=~1, share=TRUE)
  p.c.Time.share=list(formula=~Time, share=TRUE)
  
  # Candidate Models
  mod_behavior<-mark(data.processed,data.ddl,model.parameters=list(c=c.dot, p=p.dot,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.c.share<-mark(data.processed,data.ddl,model.parameters=list(p=p.c.share,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.c.Time.share<-mark(data.processed,data.ddl,model.parameters=list(p=p.c.Time.share,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.Time<-mark(data.processed,data.ddl,model.parameters=list(c=c.dot, p=p.Time,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.Time<-mark(data.processed,data.ddl,model.parameters=list(c=c.Time, p=p.dot,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.c.Time<-mark(data.processed,data.ddl,model.parameters=list(c=c.Time, p=p.Time,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  
  return(collect.models())
  
}
behavior.results=behavior.models(data.processed)

(behavior.table<-kable(model.table(model.list = behavior.results, sort = TRUE, adjust = TRUE,
                                   pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                   model.name = FALSE)))

# Step 4: Recapture -------------------------------------------------------
# Step 4
# Fit recapture models
recapture.models=function(data.processed)
{
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  #
  # Define parameter models
  #
  c.dot=list(formula=~1)
  c.Time=list(formula=~Time)
  c.session=list(formula=~session)
  c.session.Area=list(formula=~session+AREA)
  # c.Prebait=list(formula=~Prebait)
  c.Area=list(formula=~AREA)
  # c.Area.Prebait=list(formula=~AREA+Prebait)
  c.Time.Area=list(formula=~Time+AREA)
  p.global=list(formula=~GRID*Time)
  S.global=list(formula=~GRID*time)
  Gamma.Prime.best=list(formula=~1,fixed=0)
  Gamma.Double.Prime.best=list(formula=~1,fixed=0)
  # Candidate Models
  mod_c.dot<-mark(data.processed,data.ddl,model.parameters=list(c=c.dot, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.Time<-mark(data.processed,data.ddl,model.parameters=list(c=c.Time, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.session<-mark(data.processed,data.ddl,model.parameters=list(c=c.session, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.session.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.session.Area, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  # mod_c.Prebait<-mark(data.processed,data.ddl,model.parameters=list(c=c.Prebait, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.Area, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  # mod_c.Area.Prebait<-mark(data.processed,data.ddl,model.parameters=list(c=c.Area.Prebait, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_c.Time.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.Time.Area, p=p.global,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  
  
  return(collect.models())
  
}
recapture.results=recapture.models(data.processed)

(recapture.table<-kable(model.table(model.list =recapture.results, sort = TRUE, adjust = TRUE,
                                    pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                    model.name = FALSE)))

# Step 5: Capture Models --------------------------------------------------
# Step 5
# Fit capture models splitsession
capture.models=function(data.processed)
{
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  
  #
  # Define parameter models
  #
  
  c.best=list(formula=~Time+AREA)
  S.global=list(formula=~GRID*time)
  Gamma.Prime.best=list(formula=~1,fixed=0)
  Gamma.Double.Prime.best=list(formula=~1,fixed=0)
  
  p.dot=list(formula=~1)
  p.Time=list(formula=~Time)
  p.session=list(formula=~session)
  p.session.Area=list(formula=~session+AREA)
  p.Area=list(formula=~AREA)
  p.Time.Area=list(formula=~Time+AREA)
  
  # Candidate Models
  mod_p.dot<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.dot,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.Time<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.Time,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.session<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.session,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.session.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.session.Area,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.Area,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_p.Time.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.Time.Area,S=S.global,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  
  
  return(collect.models())
  
}
capture.results=capture.models(data.processed)

(capture.table<-kable(model.table(model.list = capture.results, sort = TRUE, adjust = TRUE,
                                  pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                  model.name = FALSE)))

# Step 6: Survival --------------------------------------------------------


#Step 6
#Fit Survival models
Survival.models=function(data.processed)
{
  data.ddl=make.design.data(data.processed)
  data.ddl<-fix.params(data.ddl)
  #
  # Define parameter models
  #
  c.best=list(formula=~Time+AREA)
  p.best=list(formula=~session+AREA)
  Gamma.Prime.best=list(formula=~1,fixed=0)
  Gamma.Double.Prime.best=list(formula=~1,fixed=0)
  
  S.dot=list(formula=~1)
  S.session=list(formula=~time)
  S.Area=list(formula=~AREA)
  S.session.Area=list(formula=~time+AREA)
  
  
  # Candidate Models
  mod_S.dot<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.best,S=S.dot,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_S.session<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.best,S=S.session,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_S.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.best,S=S.Area,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  mod_S.session.Area<-mark(data.processed,data.ddl,model.parameters=list(c=c.best, p=p.best,S=S.session.Area,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
  
  
  return(collect.models())
  
}
Survival.results=Survival.models(data.processed)

(Survival.table<-kable(model.table(model.list = Survival.results, sort = TRUE, adjust = TRUE,
                                   pf = 1, use.lnl = FALSE, use.AIC = FALSE,
                                   model.name = FALSE)))

# Stage 2 of the Sequential Modelling Process ------------------------------


dater1<-import.chdata("C:/Users/weldy/Projects/T_WILSON_SMALL_MAMMALS/CAPTURE_EVALUATION/data/inp.txt",header=TRUE, use.comments = TRUE)
dater1<-subset(dater1,AREA!=1)
dater2<-import.chdata("C:/Users/weldy/Projects/T_WILSON_SMALL_MAMMALS/CAPTURE_EVALUATION/data/8inp.txt",header=TRUE, use.comments = TRUE)
dater2<-subset(dater2,AREA!=1)
dater3<-import.chdata("C:/Users/weldy/Projects/T_WILSON_SMALL_MAMMALS/CAPTURE_EVALUATION/data/4inp.txt",header=TRUE, use.comments = TRUE)
dater3<-subset(dater3,AREA!=1)

dater1<-subset(dater1,ch!="000000000000000000000000000000000000000000000000000000000000000000000000")
dater2<-subset(dater2,ch!="000000000000000000000000000000000000000000000000")
dater3<-subset(dater3,ch!="000000000000000000000000")

interval1<-c(rep(c(0,0,0,0,0,0,0,0,0,0,0,1),5),c(0,0,0,0,0,0,0,0,0,0,0))
interval2<-c(rep(c(0,0,0,0,0,0,0,1),5),c(0,0,0,0,0,0,0))
interval3<-c(rep(c(0,0,0,1),5),c(0,0,0))

data.processed1=process.data(dater1,model="RDHuggins",time.intervals=interval1,groups=c("AREA","GRID"))
data.processed2=process.data(dater2,model="RDHuggins",time.intervals=interval2,groups=c("AREA","GRID"))
data.processed3=process.data(dater3,model="RDHuggins",time.intervals=interval3,groups=c("AREA","GRID"))

data.processed1$group.covariates$AREA<-relevel(data.processed1$group.covariates$AREA,ref='3')
data.processed2$group.covariates$AREA<-relevel(data.processed2$group.covariates$AREA,ref='3')
data.processed3$group.covariates$AREA<-relevel(data.processed3$group.covariates$AREA,ref='3')




data.ddl1=make.design.data(data.processed1)
data.ddl1<-fix.params(data.ddl1)
data.ddl2=make.design.data(data.processed2)
data.ddl2<-fix.params(data.ddl2)
data.ddl3=make.design.data(data.processed3)
data.ddl3<-fix.params(data.ddl3)

#
# Define parameter models
#

p.best=list(formula=~session+AREA)
S.best=list(formula=~time+AREA)
c.best=list(formula=~Time+AREA)
Gamma.Prime.best=list(formula=~1,fixed=0)
Gamma.Double.Prime.best=list(formula=~1,fixed=0)

# Run Candidate Models
mod_best.1<-mark(data.processed1,data.ddl1,model.parameters=list(c=c.best, p=p.best,S=S.best,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
mod_best.2<-mark(data.processed2,data.ddl2,model.parameters=list(c=c.best, p=p.best,S=S.best,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))
mod_best.3<-mark(data.processed3,data.ddl3,model.parameters=list(c=c.best, p=p.best,S=S.best,GammaDoublePrime=Gamma.Double.Prime.best,GammaPrime=Gamma.Prime.best))