Mark Recapture Distance Sampling, estimating trackline detection probability p(0) for additional variables such as Beaufort Sea States

Question

Shipboard line-transect distance sampling methods are commonly used to estimate populations of cetaceans. One of the assumptions of these methods is that all animals on the trackline are detected, but this is not always the case (especially for cryptic species). Mark Recapture Distance Sampling (MRDS) is a method that uses a secondary observation team to identify the fraction of animals detected on the trackline when it is suspected that animals have been missed. During typical visual line transect methods, MRDS use a secondary visual observation team to understand the fraction of animals on the trackline that are missed by the primary team. In this paper (https://link.springer.com/article/10.1007/s10651-020-00443-7), acoustic detections from a towed hydrophone array are used as the secondary 'observation' platform.

One component of the methods presented in this paper estimates the trackline detection probability p(0) for different levels of an additional variable such as Beaufort Sea States (as presented in Table 2). This was not part of the MRDS output in R. How can these methods be replicated? Specifically, if I have completed MRDS analysis for my dataset, how can I estimate trackline detection probability for individual Beaufort sea states (or other variables with multiple levels)?

Answer 1

Great you are using our methods. I have extracted the code I used for making the p0 predictions in Tables 2 and 4 for our paper. I hope this helps. Please let me know how you get on with this.

Cheers, Cornelia

### MRDS visual vs acoustic
# code written by Cornelia Oedekoven, including modified functions from mrds package

##########################################################
library(Distance)
library(mrds)

##########################################################
# Predicting p(0) in the trial configuration
# p(0) estimates per beaufort state
newdata.trial<-data.frame(distance=rep(0,6),beaufort=c(0:5))
# predict will give you the average detection probability p
# not p(0)
# doesn't matter which distance you use
g0.trial<-round(predict(object = best.trial,newdata = newdata.trial,integrate = F)$fitted, 3)
g0.trial

# average p from relative det fct g(y)
summary(best.trial)$ds.summary$average.p

# p(0) for the different beaufort states 
# => to obtain p(0) for the different beaufort states we take the p predicitons per beaufort 
# and divide them by the p from the relative detection function
p0.beaufort<-g0.trial/summary(best.trial)$ds.summary$average.p
# Table 2 (Rankin et al. 2020)
p0.beaufort

###################################################################
#### Predicting p0 in the independent observer configuration
#' @param object ddf.io object
#' @param newdata New data for making predictions
my.predict.io <-
function(object,newdata=NULL,compute=FALSE,int.range=NULL,...){

  model <- object
  if(is.null(newdata)){
    xmat <- model$mr$mr$data
  }else{
    compute <- TRUE
    xmat <- newdata
  }
  
  xmat$distance <- 0
  ddfobj <- model$ds$ds$aux$ddfobj
  
  # for gamma models need to find where p(x)=1 (apex), set that as distance
  if(ddfobj$type=="gamma"){
    xmat$distance <- rep(apex.gamma(ddfobj),2)
  }
  
  # calculate ps for each part of the model
  xmat$offsetvalue <- 0
  p.0 <- predict(model$mr,newdata=xmat,integrate=FALSE,compute=compute)
  if(is.null(newdata)){
    pdot <- predict(model$ds,esw=FALSE,compute=compute,
                    int.range=int.range)$fitted
  }else{
    pdot <- predict(model$ds,newdata=newdata[newdata$observer==1,],
                    esw=FALSE,compute=compute,int.range=int.range)$fitted
  }
  
  return(list(p.0 = p.0, pdot = pdot))
}

# example with Beaufort (Table 4 from Rankin et al. 2020) 
bft=c(0:5)
newdata.io<-data.frame(distance=rep(0,12),beaufort=rep(bft,each=2),observer=rep(c(1,2),times=6))
io.pred<-data.frame(matrix(NA,3,(length(bft)+1)))
colnames(io.pred)<-c("obs",bft)
io.pred[,1]<-c("obs1","obs2","both")
io.pred[1,2:(length(bft)+1)]<-my.predict.io(best.io,newdata.io)$p.0$p1
io.pred[2,2:(length(bft)+1)]<-my.predict.io(best.io,newdata.io)$p.0$p2
io.pred[3,2:(length(bft)+1)]<-my.predict.io(best.io,newdata.io)$p.0$fitted
print(io.pred, row.names = F)