Data=read.csv('QMG 2014 Detections Distance Double-observer.csv') Delta_measure=-(Data$DISTFRONT[which(Data$groupFLY==0)]-Data$DISTBACK[which(Data$groupFLY==0)]) Delta_move = -(Data$DISTFRONT[which(Data$groupFLY==1)]-Data$DISTBACK[which(Data$groupFLY==1)]) library(ggplot2) Plot_df = data.frame(Type = c(rep("Flying",length(Delta_move)),rep("Not flying",length(Delta_measure))), Error = c(Delta_move,Delta_measure)) Dist_error_hist = ggplot(Plot_df)+geom_histogram(aes(Error))+facet_grid(.~Type)+ylab('Number of waterfowl groups')+ theme(axis.text=element_text(size=16),axis.title=element_text(size=16),strip.text.x = element_text(size = 16)) png("Dist_error_hists.png") Dist_error_hist dev.off() Data$DISTFRONT[Data$DISTFRONT==6]=NA Data$Dist = Data$DISTFRONT Data$Dist[is.na(Data$Dist)]=Data$DISTBACK[is.na(Data$Dist)] Data$Dist[Data$Dist==6]=NA # So, only include observations for which distance of first observer was in bins 1-5 OR first observer missed it but second observer saw it in bins 1-5; i.e. observer 2 distance bin 6 is still included as a detection of first observer saw it in 1-5 Data=Data[-which(is.na(Data$Dist)),] Data$Species=as.character(Data$Species) Data$Species[Data$Species=="UK LOON"]="UKLOON" Data$Species = factor(Data$Species) Which.sp = unique(Data$Species)[which(tabulate(Data$Species)>=20)] #which species seen 20 or more times Data=Data[Data$Species %in% Which.sp,] Data$Species = factor(as.character(Data$Species)) Species.list = unique(Data$Species) n.species=length(Species.list) Species=Data$Species Data$Species=as.character(Data$Species) for(isp in 1:n.species)Data$Species[Data$Species==Species.list[isp]]=isp Data$Species=as.numeric(Data$Species) #Data$Species = factor(as.character(Data$Species)) #reset species levels for factor var n.transects = 1 n.obs = nrow(Data) ch_split <- function(x){ as.numeric(c(substr(x,1,1),substr(x,2,2))) } Data$CH = as.character(Data$CH) Data$CH[Data$CH=="1"]="01" Obs = as.vector(sapply(Data$CH,"ch_split")) Dat = data.frame(Transect = rep(1,n.obs*2),Match = rep(c(1:n.obs),each=2),Observer=rep(c(1,2),n.obs),Obs=Obs,Species=rep(Data$Species,each=2),Distance = rep(Data$Dist,each=2),Group=rep(Data$GROUPSZ,each=2),Fly=rep(Data$groupFLY,each=2)) Dat_LP = Dat Dat_LP[,"Distance"]=1 #conduct double observer analysis using L-P and stratified L-P w Huggins in RMARK- note #estimates are only for # of groups; estimates of # of individuals would need to be made using a #bootstrap with group/cluster size in the numerator of the H-T procedure Dat_MARK = Data Dat_MARK$ch = Dat_MARK$CH Dat_MARK$Cluster = Dat_MARK$GROUPSZ Dat_MARK$Species = Species Dat_MARK$NewSpecies Dat_MARK$Obs11 = 0 Dat_MARK$Obs12 = as.numeric(substr(Dat_MARK$CH,1,1)) Dat_MARK$Dist21 = (Dat_MARK$Dist==2)*1 Dat_MARK$Dist31 = (Dat_MARK$Dist==3)*1 Dat_MARK$Dist41 = (Dat_MARK$Dist==4)*1 Dat_MARK$Dist51 = (Dat_MARK$Dist==5)*1 Dat_MARK$Dist22 = (Dat_MARK$Dist==2)*1 Dat_MARK$Dist32 = (Dat_MARK$Dist==3)*1 Dat_MARK$Dist42 = (Dat_MARK$Dist==4)*1 Dat_MARK$Dist52 = (Dat_MARK$Dist==5)*1 Dat_MARK$FlyDist21 = (Dat_MARK$Dist==2)*Dat_MARK$groupFLY Dat_MARK$FlyDist31 = (Dat_MARK$Dist==3)*Dat_MARK$groupFLY Dat_MARK$FlyDist41 = (Dat_MARK$Dist==4)*Dat_MARK$groupFLY Dat_MARK$FlyDist51 = (Dat_MARK$Dist==5)*Dat_MARK$groupFLY Dat_MARK$FlyDist22 = (Dat_MARK$Dist==2)*Dat_MARK$groupFLY Dat_MARK$FlyDist32 = (Dat_MARK$Dist==3)*Dat_MARK$groupFLY Dat_MARK$FlyDist42 = (Dat_MARK$Dist==4)*Dat_MARK$groupFLY Dat_MARK$FlyDist52 = (Dat_MARK$Dist==5)*Dat_MARK$groupFLY Dat_MARK$Ddist = abs(Dat_MARK$DISTBACK-2) Dat_MARK$Ddist[is.na(Dat_MARK$Ddist)]=0 library('RMark') HA_proc = process.data(Dat_MARK,model="Huggins",groups="Species") HA_ddl = make.design.data(HA_proc) p_model = list(formula=~1,share=TRUE) HA_mod_dot = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Cluster,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~groupFLY,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~time,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster+time,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Cluster+time,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+time,share=TRUE) HA_mod_SOG = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Cluster+time+groupFLY,share=TRUE) HA_mod_SOGF = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster+time+groupFLY,share=TRUE) HA_mod_SOGF = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster+time+groupFLY+Dist2:time+Dist3:time+Dist4:time+Dist5:time,share=TRUE) HA_mod_SOGFD = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster+time+groupFLY+Dist2:time+Dist3:time+Dist4:time+Dist5:time+FlyDist2:time+FlyDist3:time+FlyDist4:time+FlyDist5:time,share=TRUE) HA_mod_SOGFDint = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Species+Cluster+time+Dist2:time+Dist3:time+Dist4:time+Dist5:time,share=TRUE) HA_mod_SOGFD = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model)) p_model = list(formula=~Cluster+time+groupFLY+Dist2:time+Dist3:time+Dist4:time+Dist5:time,share=TRUE) HA_mod_SOGFD = mark(HA_proc,HA_ddl,model.parameters=list(p=p_model))