############################################################################################
#                                                                                          #
#                                        SOFT REQUIRED                                     #
#                                                                                          #
############################################################################################

# mark program should be installed on the system for this script to work
# RMark package should be installed on the system

############################################################################################
#                                                                                          #
#                                          PACKAGES                                        #
#                                                                                          #
############################################################################################

### packages used in analyses

library(ggplot2)
library(RMark)

############################################################################################
#                                                                                          #
#                                    NEST DAILY SURVIVAL RATE                              #
#                                                                                          #
############################################################################################

# removing all variables
rm(list = ls(all = TRUE))

# reading data table; has to be in the working directory

lsrk <- read.table("LSvsRK.txt", header=T, sep=";",na.strings='')
lsrk

# adding year as factor variable 

lsrk$YearF = factor(lsrk$Year)

# A model of constant daily survival rate (DSR)

S.lsrk=mark(lsrk,nocc=48,model="Nest", model.parameters=list(S=list(formula=~1)))
S.lsrk$results$real

# model with species as factor

S.Sp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~Sp)), groups=c("Sp"))
S.Sp$results$real

# model with year as factor

S.YearF=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF)), groups=c("YearF"))
S.YearF$results$real

# model with elevation as covariate

S.ElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~ElevMEAN)))
S.ElevMEAN$results$real

# model with evelation as factor

S.Elev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~Elev60)), groups=c("Elev60"))
S.Elev60$results$real



################### selection of the best model #################################

run.lsrkF=function()
{

# A model of constant daily survival rate (DSR)
   S.lsrk =mark(lsrk,model="Nest",nocc=48)

# model with year as factor

   S.YearF=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF)), groups=c("YearF"))

# model with elevation as covariate

  S.ElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~ElevMEAN)))

# model with evelation as factor

  S.Elev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~Elev60)), groups=c("Elev60"))

# model with species as factor

  S.Sp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~Sp)), groups=c("Sp"))

# model with year and species as factors

   S.YearF.Sp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Sp)), groups=c("YearF","Sp"))

# model with year and species as factors, elevation as covariate

   S.YearF.Sp.ElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Sp+ElevMEAN)), groups=c("YearF","Sp"))

# model with year as factor and elevation as covariate

   S.YearF.ElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+ElevMEAN)), groups=c("YearF"))

# model with year, species and evelation as factors

   S.YearF.Sp.Elev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Sp+Elev60)), groups=c("YearF","Sp","Elev60"))

# model with year and elevation as factors

   S.YearF.Elev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Elev60)), groups=c("YearF","Elev60"))

# models including interactions

# Year * species

   S.YearFxSp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Sp)), groups=c("YearF","Sp"))

# Year * species * evelation

   S.YearFxSpxElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Sp*ElevMEAN)), groups=c("YearF","Sp"))

# Year * species + evelation

   S.YearFxSp.ElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Sp+ElevMEAN)), groups=c("YearF","Sp"))

# Year + species * elevation

   S.YearF.SpxElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Sp*ElevMEAN)), groups=c("YearF","Sp"))

# Year * elevation + species

   S.YearFxElevMEAN.Sp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*ElevMEAN+Sp)), groups=c("YearF","Sp"))

# Year * elevation

   S.YearFxElevMEAN=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*ElevMEAN)), groups=c("YearF"))


# Year * species * âûñîòà, all as factors

   S.YearFxSpxElev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Sp*Elev60)), groups=c("YearF","Sp","Elev60"))

# Year * species + elevation, all as factors

   S.YearFxSp.Elev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Sp+Elev60)), groups=c("YearF","Sp","Elev60"))

# Year + species * elevation, all as factors

   S.YearF.SpxElev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF+Sp*Elev60)), groups=c("YearF","Sp","Elev60"))

# Year * elevation + species, all as factors

   S.YearFxElev60.Sp=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Elev60+Sp)), groups=c("YearF","Sp","Elev60"))

# Year * elevation, all as factors

   S.YearFxElev60=mark(lsrk,model="Nest",nocc=48, model.parameters=list(S=list(formula=~YearF*Elev60)), groups=c("YearF","Elev60"))

   return(collect.models())
}

# run defined models

lsrkF.results=run.lsrkF()

# displaying table of model results

lsrkF.results


# displaying details for the best 3 models

lsrkF.results$S.YearF.Sp			# print MARK output to designated text editor
lsrkF.results$S.YearF.Sp$results$beta	# view estimated beta’s in R
lsrkF.results$S.YearF.Sp$results$real	# view estimated DSR estimate in R

lsrkF.results$S.YearFxSp			# print MARK output to designated text editor
lsrkF.results$S.YearFxSp$results$beta	# view estimated beta’s in R
lsrkF.results$S.YearFxSp$results$real	# view estimated DSR estimate in R

lsrkF.results$S.YearF.Sp.ElevMEAN			# print MARK output to designated text editor
lsrkF.results$S.YearF.Sp.ElevMEAN$results$beta	# view estimated beta’s in R
lsrkF.results$S.YearF.Sp.ElevMEAN$results$real	# view estimated DSR estimate in R

### daily survival rate graph

survival_daily<-read.csv("Real_Function_Parameters_of_the_best_survival_model.csv")
str(survival_daily)
ggplot(data=survival_daily,aes(y=estimate,x=as.factor(year),fill=sp))+
   geom_point(shape=21,size=5)+
   geom_line()+
   scale_y_continuous(limits=c(0.9,1))+
   theme_classic()

# cleanup

cleanup(ask=FALSE)
list.files()