##This is the simpler analysis using standard double observer analysis, no covariates. ##remove obsservations before 0.2 meters B <- 1e4 #size of bootstrap, increase for more precision. manuscript uses 10^4 library(dplyr) library(MASS) library(stringr) library(ggplot2) library(rgdal) library(Distance) library(spatstat) library(mgcv) library(RColorBrewer) library(kableExtra) library(numDeriv) library(nloptr) library(dsm) source("ZebraFuncs_Vanilla.R") #numerous functions used to format the raw data. ###This chunk reads in the datasets, the formats them for our analyses### # read in the datasets Hab <- read.csv(file="HabitatCharacterizat_DATA_2017-10-12_0821.csv", na.string="9999") #habitat data Trans <- read.csv(file="TransectData_DATA_2017-10-23_1150.csv", na.string="9999") #transect data ZM <- read.csv(file="ZebraMusselEncouter_DATA_2017-10-12_0823.csv", na.string="9999") #encounter data #renames some transect columns for consistency with encounter and habitat datasets. Finally subsets data to include only lake Sylvia and lake Burgan Trans <- Trans %>% rename(Region.Label = lake_name) %>% rename(gps_transect_n="gps_location_n") %>% rename(gps_transect_w="gps_location_w") %>% subset(Region.Label == "Sylvia" | Region.Label == "Burgen") #renames encounter entries to be consistent with what is needed for distance analyses, create a record_id for each detection, converts detection distance to meters. ZM <- ZM %>% rename(distance=perp_distance) %>% rename(encounter_id=record_id) %>% rename(size=number_mussels) %>% mutate(record_id = as.factor(str_split(encounter_id, pattern="_E", n=2, simplify = TRUE)[,1])) %>% mutate(distance=distance/100) #impute missing values #Trans[which(is.na(Trans$clarity)),]$clarity <- mean(Trans$clarity, na.rm=TRUE) # interpolates missing Trans data in lake sylvia Trans <- interpolate(Trans) #' merges transect information and encounter information ZM <- merge(Trans, ZM, by='record_id') #' Determine which encounters are associated with the second observers. Link those encounters to observers 5 & 6. for(i in 1:length(ZM$encounter_id)) { temp <- str_sub(ZM$encounter_id[i], start=nchar(as.character(ZM$encounter_id[i]))-1, end=nchar(as.character(ZM$encounter_id[i]))) if(temp == "_2") { ZM$observer___1[i] <- 0 ZM$observer___2[i] <- 0 ZM$observer___5[i] <- 1 ZM$observer___6[i] <- 1 } } #this ection links habitat type to each encounter. ZM <- cbind(ZM, substrate1=rep(NA, dim(ZM)[1]), substrate2=rep(NA, dim(ZM)[1]), substrate3=rep(NA, dim(ZM)[1]), substrate4=rep(NA, dim(ZM)[1]), substrate5=rep(NA, dim(ZM)[1]), substrate6=rep(NA, dim(ZM)[1]), substrate7=rep(NA, dim(ZM)[1]), plantcover=rep(NA, dim(ZM)[1]), Effort=rep(NA, dim(ZM)[1]), Area=rep(NA, dim(ZM)[1])) #fix one mislabeled entry Hab$record_id <- recode(Hab$record_id, "21004900_T03_30"="21004900_T03_D30") # cleanup hab dataset such as adding in effort for transects without any observations temp <- Link.HabCounts(Hab, ZM) HabZM <- temp[[1]] #linked hab/encounter data Hab <- temp[[2]] #original habitat data temp <- matrix(NA, dim(HabZM)[1], 2) for(i in 1:dim(HabZM)[1]) { temp[i,] <- detection.location(HabZM[i,]) } HabZM <- mutate(HabZM, gps_detect_w=temp[,1], gps_detect_n=temp[,2]) #encounter data for each lake SylviaDat <- filter(HabZM, Region.Label == "Sylvia") BurgenDat <- filter(HabZM, Region.Label == "Burgen") BurgenDat <- mutate(BurgenDat, observer=(BurgenDat$observer___1==1) + 2*(BurgenDat$observer___5==1)) #create observer variable for Lake Burgen #create segments out of transect data. Used in count analysis. SegDat <- create.Seg(Hab, Trans, lakename="Burgen") #match observations made by each observer, using the distance criterion described in manuscript. BurgenDDF <- match.function(BurgenDat, horDis=0.2) ObsDat <- create.Obs(BurgenDat, SegDat=SegDat, lakename="Burgen") #makes the count data in each segment #mlv(x=BurgenDDF$distance, method="lientz", bw=0.05) dthresh <- 0.2 BurgenDDF <- BurgenDDF[-which(BurgenDDF$distance <= dthresh),] BurgenDDF$distance <- BurgenDDF$distance - dthresh #rescale for distance model SylviaDat <- SylviaDat[-which(SylviaDat$distance <= dthresh),] SylviaDat$distance <- SylviaDat$distance - dthresh ##Sylvia detection analysis ############### varS <- var(SylviaDat$size)/length(SylviaDat$size) eS <- mean(SylviaDat$size) SylviaDis <- data.frame(distance=SylviaDat$distance, object=1:length(SylviaDat$distance), observer=1, detected=1) SylviaDis.model <- ddf(method="ds", dsmodel=~cds(key="hn", formula=~1), data=SylviaDis, meta.data=list(width=1-dthresh)) Sylvia.phat <- summary(SylviaDis.model)$average.p Sylvia.phat.se <- summary(SylviaDis.model)$average.p.se[1,1] ###################################### ###now estimate the detection function for lake Burgan BurganUnique <- BurgenDDF Burgan.Trans <- Trans[which(Trans$Region.Label == "Burgen"),] for(i in 1:64) { ind <- which(BurganUnique == i) if(length(ind) > 1) { BurganUnique[-ind[-1],] } } #fit the sight-resight component of the detection function ddf.dsVanilla <- ddf(method="io", dsmodel=~cds(key="hn", formula=~ 1), mrmodel=~glm(link="logit", formula=~1), data=BurgenDDF, meta.data=list(width=1-dthresh)) Burgan.phat <- summary(ddf.dsVanilla)$average.p Burgan.phat.se <- summary(ddf.dsVanilla)$average.p.se[1,1] count.vec <- vector('numeric', max(Burgan.Trans$transect_no)) for(i in 1:length(count.vec)) { count.vec[i] <- length(unique(BurganUnique[which(BurganUnique$transect_no == i),]$object)) } area.vec <- vector('numeric', max(Burgan.Trans$transect_no)) for(i in 1:length(area.vec)) { area.vec[i] <- sum(SegDat[SegDat$transect_no==i,]$Effort)*1.6 } Burgan.dhat <- sum(count.vec, na.rm=T)/sum(area.vec)/Burgan.phat K <- length(count.vec) l <- area.vec/1.6 L <- sum(l) n <- count.vec N <- sum(count.vec) var.n <- L/(K-1)*sum(l*(n/l - N/L)^2) Burgan.dhat.se.phat <- sqrt(Burgan.dhat^2*Burgan.phat.se^2/Burgan.phat^2) Burgan.dhat.se <- sqrt(Burgan.dhat^2*(var.n/N^2 + Burgan.phat.se^2/Burgan.phat^2)) save(list=c('SylviaDat', 'Sylvia.phat', 'Sylvia.phat.se', 'SylviaDis.model', 'Burgan.phat', 'Burgan.phat.se', 'Burgan.dhat', 'Burgan.dhat.se', 'Burgan.dhat.se.phat', 'ddf.dsVanilla'), file="ZeebDensityResults_Vanilla.Rdata")