############################## # # # 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))