##ARCTIC PEREGRINE ANALYSIS -- cliff abundance using unmarked in R 

#Set working directory to folder containing the CSV files
#setwd("C:/Users/Username/Desktop/folder_containing_csvs/")


##abundance and covariate input files
obsmatrix=read.csv("obsmatrix.csv",header=T)
visitmatrix=read.csv("visitmatrix.csv",header=T)
sitecovs=read.csv("sitecovs.csv",header=T)
productivitymatrix=read.csv("productivitymatrix.csv",header=T)
precipmatrix=read.csv("precipmatrix.csv",header=T)
meltdatematrix=read.csv("meltdatematrix.csv",header=T)
tmaxmatrix=read.csv("tmaxmatrix.csv",header=T)
abundancematrix=read.csv("abundancematrix.csv",header=T)
yearsmatrix=read.csv("yearsmatrix.csv",header=T)
yearslogmatrix=read.csv("yearslogmatrix.csv",header=T)
yearsthresholdmatrix=read.csv("yearsthresholdmatrix.csv",header=T)
gyrdistancematrix=read.csv("gyrdistancematrix.csv",header=T)
pdomatrix=read.csv("pdomatrix.csv",header=T)

####covariate years as a factor
years = as.character(1982:2002)
years = matrix(years, nrow=74, 21, byrow=TRUE)

##Install package "unmarked" if you haven't yet
#install.packages('unmarked')

####call to package unmarked in R; modifying input files to remove first column
library(unmarked)
ymod=obsmatrix[,2:43]
visitmod=visitmatrix[,2:43]
productivitymod=productivitymatrix[,2:22]
precipmod=precipmatrix[,2:22]
meltdatemod=meltdatematrix[,2:22]
tmaxmod = tmaxmatrix[,2:22]
abundancemod = abundancematrix[,2:22]
yearsmod = yearsmatrix[,2:22]
yearslogmod = yearslogmatrix[,2:22]
yearsthresholdmod = yearsthresholdmatrix[,2:22]
gyrdistancemod = gyrdistancematrix[,2:22]
pdomod = pdomatrix[,2:22]

####creation of dataframe in unmarked
nestabundanceinput = unmarkedFramePCO(y=ymod, siteCovs=sitecovs, obsCovs=list(visit=visitmod), yearlySiteCovs=list(year=years, yearlinear=yearsmod, yearthreshold=yearsthresholdmod, yearlog=yearslogmod, productivityt1=productivitymod, precip=precipmod, meltdate=meltdatemod, tmax=tmaxmod, abundancet1=abundancemod, gyrdistance=gyrdistancemod, pdo=pdomod), numPrimary=21)
summary(nestabundanceinput)

####final model list to be evaluated
m1 = pcountOpen(lambdaformula=~1, gammaformula=~1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m1)
m2 = pcountOpen(lambdaformula=~geology, gammaformula=~aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m2)
m3 = pcountOpen(lambdaformula=~height, gammaformula=~aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m3)
m4 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m4)
m5 = pcountOpen(lambdaformula=~geology, gammaformula=~geology + aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m5)
m6 = pcountOpen(lambdaformula=~height, gammaformula=~geology + aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m6)
m7 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~geology + aspectcat, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m7)
m8 = pcountOpen(lambdaformula=~geology, gammaformula=~meltdate+pdo, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m8)
m9 = pcountOpen(lambdaformula=~height, gammaformula=~meltdate+pdo, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m9)
m10 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~meltdate+pdo, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m10)
m11 = pcountOpen(lambdaformula=~geology, gammaformula=~pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m11)
m12 = pcountOpen(lambdaformula=~height, gammaformula=~pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m12)
m13 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m13)
m14 = pcountOpen(lambdaformula=~geology, gammaformula=~abundancet1+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m14)
m15 = pcountOpen(lambdaformula=~height, gammaformula=~abundancet1+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m15)
m16 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~abundancet1+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m16)
m17 = pcountOpen(lambdaformula=~geology, gammaformula=~gyrdistance+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m17)
m18 = pcountOpen(lambdaformula=~height, gammaformula=~gyrdistance+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m18)
m19 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~gyrdistance+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m19)
m20 = pcountOpen(lambdaformula=~geology, gammaformula=~meltdate+pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m20)
m21 = pcountOpen(lambdaformula=~height, gammaformula=~meltdate+pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m21)
m22 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~meltdate+pdo+productivityt1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m22)
m23 = pcountOpen(lambdaformula=~geology, gammaformula=~meltdate+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m23)
m24 = pcountOpen(lambdaformula=~height, gammaformula=~meltdate+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m24)
m25 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~meltdate+pdo+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m25)
m26 = pcountOpen(lambdaformula=~geology, gammaformula=~pdo+productivityt1+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m26)
m27 = pcountOpen(lambdaformula=~height, gammaformula=~pdo+productivityt1+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m27)
m28 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~pdo+productivityt1+tmax, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m28)
m29 = pcountOpen(lambdaformula=~geology, gammaformula=~1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m29)
m30 = pcountOpen(lambdaformula=~height, gammaformula=~1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m30)
m31 = pcountOpen(lambdaformula=~height + waterarea, gammaformula=~1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m31)
m32 = pcountOpen(lambdaformula=~1, gammaformula=~aspectcat, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m32)
m33 = pcountOpen(lambdaformula=~1, gammaformula=~geology + aspectcat, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m33)
m34 = pcountOpen(lambdaformula=~1, gammaformula=~meltdate + pdo, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m34)
m35 = pcountOpen(lambdaformula=~1, gammaformula=~pdo + productivityt1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m35)
m36 = pcountOpen(lambdaformula=~1, gammaformula=~abundancet1+pdo+tmax, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m36)
m37 = pcountOpen(lambdaformula=~1, gammaformula=~gyrdistance+pdo+tmax, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m37)
m38 = pcountOpen(lambdaformula=~1, gammaformula=~meltdate+pdo+productivityt1, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m38)
m39 = pcountOpen(lambdaformula=~1, gammaformula=~meltdate+pdo+tmax, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m39)
m40 = pcountOpen(lambdaformula=~1, gammaformula=~pdo+productivityt1+tmax, omegaformula=~1, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m40)
m41 = pcountOpen(lambdaformula=~1, gammaformula=~1, omegaformula=~abundancet1+pdo+tmax, pformula=~1, mixture="ZIP", data=nestabundanceinput, K=20, dynamics="constant")
summary(m41)
models = fitList('m1'=m1,'m2'=m2,'m3'=m3,'m4'=m4,'m5'=m5,'m6'=m6','m7'=m7,'m8'=m8,'m9'=m9,'m10'=m10,'m11'=m11,'m12'=m12,'m13'=m13,'m14'=m14,'m15'=m15,'m16'=m16,'m17'=m17,'m18'=m18,'m19'=m19,'m20'=m20,'m21'=m21,'m22'=m22,'m23'=m23,'m24'=m24,'m25'=m25,'m26'=m26,'m27'=m27,'m28'=m28,'m29'=m29,'m30'=m30,'m31'=m31,'m32'=m32,'m33'=m33,'m34'=m34,'m35'=m35,'m36'=m36,'m37'=m37,'m38'=m38,'m39'=m39,'m40'=m40,'m41'=m41)
ms = modSel(models)
ms