library(unmarked)
library(AICcmodavg)

setwd("C:\\Users\\obses\\Desktop")
data1<-read.csv("colext_data_practice_1.csv", header=TRUE)

M <- (nrow(data1)) # Number of sites = m
J <- 2 # Number secondary sample periods
T <- 3 # Number primary sample periods

# matrix(data, nrow, ncol)

yy <- matrix(c(data1$det1.1, data1$det1.2, data1$det2.1, data1$det2.2, data1$det3.1, data1$det3.2), M, J*T) # make a matrix of the input data (detection histories across all years)
time <- matrix(c(data1$time1.1, data1$time1.2, data1$time2.1, data1$time2.2, data1$time3.1, data1$time3.2), M, J*T)
date <- matrix(c(data1$date1.1, data1$date1.2, data1$date2.1, data1$date2.2, data1$date3.1, data1$date3.2), M, J*T)
thistle<- matrix(c(data1$thistle), nrow(data1), 1, byrow=TRUE)
dist<- matrix(c(data1$dist), nrow(data1), 1, byrow=TRUE)
year <- matrix(c('1','2','3'), nrow(data1), T, byrow=TRUE)

###########

umf1 <- unmarkedMultFrame(y = yy,
                          yearlySiteCovs = list(year = year),
                          siteCovs= data.frame(thistle=thistle, dist=dist),
                          obsCovs = list(time = time, date = date), numPrimary=T)

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

p.null <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ 1, data = umf1, method="BFGS")
p.date <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ date, data = umf1, method="BFGS")
p.time <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ time, data = umf1, method="BFGS")
p.year <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ year, data = umf1, method="BFGS")

modlist1<-list(p.null=p.null, p.date=p.date, p.time=p.time, p.year=p.year)
aictab(modlist1)

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

confint(p.date, type="det")

# Graphing detection over date

nd <- data.frame(("date"=(seq(0, 1, length.out = 100))))
colnames(nd)<-c("date")
pred.p <- predict(p.date, type='det', newdata=nd)
print(predictions<-cbind(pred.p, nd))

plot(1, xlim=c(0,1), ylim=c(0.2,0.9), type="n", axes=T, xlab="Date",
     pch=20, ylab="Detection Probability", family="serif",
     cex.lab=1.25, cex.main=1.75)

lines(predictions$date, predictions$Predicted, col="black", lwd=2)
lines(predictions$date, predictions$lower, lty=2, col="black")
lines(predictions$date, predictions$upper, lty=2, col="black")

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

psi.null <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ date, data = umf1, method="BFGS")

psi.thistle <- colext(psiformula= ~thistle, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ date, data = umf1, method="BFGS")

modlist2<-list(psi.null=psi.null, psi.thistle=psi.thistle)
aictab(modlist2)

###########

confint(psi.thistle, type="psi")

##########

nd2 <- data.frame(("thistle"=(seq(0, 100, length.out = 100))))
colnames(nd2)<-c("thistle")
pred.psi <- predict(psi.thistle, type='psi', newdata=nd2)
print(predictions<-cbind(pred.psi, nd2))

plot(1, xlim=c(0,100), ylim=c(0,1), type="n", axes=T, xlab="Thistle (% cover)",
     pch=20, ylab="First-year Occupancy Probability", family="serif",
     cex.lab=1.25, cex.main=1.75)

lines(predictions$thistle, predictions$Predicted, col="black", lwd=2)
lines(predictions$thistle, predictions$lower, lty=2, col="black")
lines(predictions$thistle, predictions$upper, lty=2, col="black")

###########

gamma.null <- colext(psiformula= ~thistle, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ date, data = umf1, method="BFGS")

gamma.thistle <- colext(psiformula= ~thistle, gammaformula = ~ thistle, epsilonformula = ~ 1, pformula = ~ date, data = umf1, method="BFGS")

modlist3<-list(gamma.null=gamma.null, gamma.thistle=gamma.thistle)
aictab(modlist3)

###########

nd3 <- data.frame(("thistle"=(seq(0, 100, length.out = 100))))
colnames(nd3)<-c("thistle")
pred.gam <- predict(gamma.thistle, type='col', newdata=nd3)
print(predictions<-cbind(pred.gam, nd3))

plot(1, xlim=c(0,100), ylim=c(0,1), type="n", axes=T, xlab="Thistle (% cover)",
     pch=20, ylab="Colonization Probability", family="serif",
     cex.lab=1.25, cex.main=1.75)

lines(predictions$thistle, predictions$Predicted, col="black", lwd=2)
lines(predictions$thistle, predictions$lower, lty=2, col="black")
lines(predictions$thistle, predictions$upper, lty=2, col="black")