##############################
#                            #
#     AEC501: Ornithology    #
# Multiple Observer Sampling #
#                            #
##############################

rm(list=ls())

library(unmarked)

# Simulate independent double observer data
nSites <- 50 #number of survey points
lambda <- 10 #average number of birds per survey point 
p <- 0.3 #detection probability
cp <- c(p, p, p*p) #capture probability cell structure

set.seed(501)
N <- rpois(nSites, lambda) #Superpopulation--total across all points
PointCounts <- matrix(NA, nSites, 3)
for(i in 1:nSites) {
PointCounts[i,] <- rmultinom(1, N[i], c(cp, 1-sum(cp)))[1:3]
}

colnames(PointCounts) <- c("ObsA", "ObsB", "Both")
PointCounts <- data.frame(PointCounts)

# Take a look at our simulated data
head(PointCounts)

# Lincoln-Peterson Estimate

#PointCounts$EstN <- rep(NA, nrow(PointCounts))
#for(i in 1:nrow(PointCounts)){
#	PointCounts$EstN[i] <- (PointCounts$ObsA[i]*PointCounts$ObsB[i])/PointCounts$Both[i]
#	}

# This won't work if we have zeros in the denominator, so we'll use
# An adjusted estimate for now

# Modified Lincoln-Peterson (Chapman, 1951)
PointCounts$ModN <- rep(NA, nrow(PointCounts))
for(i in 1:nrow(PointCounts)){
	PointCounts$ModN[i] <- (((PointCounts$ObsA[i]+1)*(PointCounts$ObsB[i]+1))/
					(PointCounts$Both[i]+1))
					-1
	}

(EstN <- sum(PointCounts$ModN)) #Estimated Superpopulation
(NperPlot <- EstN/nrow(PointCounts)) #Estimated average birds per plot (lambda)

# Compare these results to our original values
sum(N)
lambda

# Let's compare these estimates to or observed values
PointCounts$ObsN <- (PointCounts$ObsA + PointCounts$ObsB) - PointCounts$Both
head(PointCounts)

### For more advanced analysis, may want to used 'unmarked' ###

rm(list=ls())
libarary(unmarked)

# Simulate independent double observer data
nSites <- 50 #number of survey points
lambda <- 10 #average number of birds per survey point 
p <- 0.3 #detection probability
cp <- c(p*(1-p), p*(1-p), p*p) #capture probability cell structure
	#Notice that this program requires a different cell structure!
	#(Observer A ONLY, Observer B ONLY, Both Observers)

set.seed(501)
N <- rpois(nSites, lambda) #Superpopulation--total across all points
PointCounts <- matrix(NA, nSites, 3)
for(i in 1:nSites) {
PointCounts[i,] <- rmultinom(1, N[i], c(cp, 1-sum(cp)))[1:3]
}

# Fit model
umf <- unmarkedFrameMPois(y=PointCounts,type="double")
fm <- multinomPois(~1 ~1, umf) #Could add covariates to this if necessary
# Estimates of fixed effects
e <- coef(fm)
exp(e[1])
plogis(e[2])
# Estimates of random effects
re <- ranef(fm, K=20)
#ltheme <- canonical.theme(color = FALSE)
#lattice.options(default.theme = ltheme)
plot(re, layout=c(10,5))