################################################################################
#
#                     WWDO Recruitment Analysis
#               2009 and 2010 Alice Capture-Recapture Data
#               Author:  Bret Collier, Texas A&M University
#               Contact:  bret@tamu.edu for questions
#
################################################################################
#Used for cleaning up workspace
rm(list=ls(all=TRUE))
#Requires RMark and MARK installed

#set working dir and read raw dataset into a appropriate frame
#setwd()
x=read.csv("Alice2009EncHist_Final_POPAN.csv", header=TRUE)
#x=read.csv("Alice2010EncHist_Final_POPAN.csv", header=TRUE)
x=data.frame(ch=paste(x$E0,x$E1,x$E2,x$E3,x$E4,x$E5,x$E6,x$E7,x$E8,x$E9,x$E10,x$E11,x$E12,x$E13, sep=""), Age=factor(x$Age), stringsAsFactors=FALSE)

#Date Windows for Enc Occasions
#1	27 Feb
#2	13 March
#3	27 March
#4	10 April
#5	24 April
#6	8 May
#7	22 May
#8	5 June
#9	19 June
#10	3 July
#11	17 July
#12	31 July
#13	14 August-End

#Load required RMark package
library(RMark)

#Process WWDO dataset for analysis
wwdo=process.data(x, model="POPAN", groups="Age")
wwdo.ddl=make.design.data(wwdo)

#Create artificial time periods grouping
wwdo.ddl=add.design.data(wwdo,wwdo.ddl,parameter="Phi",type="time",
		bins=c(0,2,4,7,11,12,13),name="Test")

#Fixing Phi Parameters for sampling periods where HY WWDO were not available in population
hy.phi1=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==1,]))
hy.phi2=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==2,]))
hy.phi3=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==3,]))
hy.phi4=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==4,]))
hy.phi5=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==5,]))
hy.phi6=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==6,]))
hy.phi7=as.numeric(row.names(wwdo.ddl$Phi[wwdo.ddl$Phi$group=="HY" & wwdo.ddl$Phi$time==7,]))
hy.phi.fix=c(hy.phi1, hy.phi2, hy.phi3, hy.phi4, hy.phi5, hy.phi6, hy.phi7)

#Fixing PENT Parameters for sampling period where HY WWDO were not available in population
hy.pent2=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==2,]))
hy.pent3=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==3,]))
hy.pent4=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==4,]))
hy.pent5=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==5,]))
hy.pent6=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==6,]))
hy.pent7=as.numeric(row.names(wwdo.ddl$pent[wwdo.ddl$pent$group=="HY" & wwdo.ddl$pent$time==7,]))
hy.pent.fix=c(hy.pent2, hy.pent3, hy.pent4, hy.pent5, hy.pent6, hy.pent7)

#####
#Real Parameter Definitions
#####
#Detection process
	p.dot=list(formula=~1)
	p.time=list(formula=~time)
	p.group=list(formula=~group)
	p.g.time=list(formula=~group:time)
	
#Survival process
	Phi.dot.fix=list(formula=~1, fixed=list(index=hy.phi.fix, value=c(0,0,0,0,0,0,0)))
	Phi.time.fix=list(formula=~time, fixed=list(index=hy.phi.fix, value=c(0,0,0,0,0,0,0)))
	Phi.age.fix=list(formula=~group, fixed=list(index=hy.phi.fix, value=c(0,0,0,0,0,0,0)))
	Phi.timeage.fix=list(formula=~time:group, fixed=list(index=hy.phi.fix, value=c(0,0,0,0,0,0,0)))
	Phi.Test=list(formula=~Test,fixed=list(index=hy.phi.fix, value=c(0,0,0,0,0,0,0)))
	
#Entry Process-always time dependent, otherwise makes no sense in my situation
	pent.time.fix=list(formula=~time, fixed=list(index=hy.pent.fix, value=c(0,0,0,0,0,0)))
	
#Model Runs
	#Note:  'N' formula should always be N=list(formula=~group)

Model.1=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.dot.fix, p=p.dot, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)
Model.2=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.time.fix, p=p.dot, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)
Model.3=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.age.fix, p=p.dot, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)
Model.4=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.timeage.fix, p=p.dot, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)
Model.5=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.timeage.fix, p=p.time, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)
Model.6=mark(wwdo, wwdo.ddl, model.parameters=list(Phi=Phi.timeage.fix, p=p.g.time, pent=pent.time.fix, N=list(formula=~group)), invisible=FALSE)

#Model selection evaluation
options(width=400)
collect.models()
Model.6