dsm_var_gam fails with strip transects

[erex]

Message to the Distance list regarding variance estimation when using strip transects with DSM. User first tried dsm_var_prop, but that failed, so tried dsm_var_gam which also failed.

Minimum reproducible example with Spotted dolphin data shipped with dsm package

library(Distance)
library(dsm)
# load the Gulf of Mexico dolphin data (see ?mexdolphins)
data(mexdolphins)
uni.mod <- ds(distdata, key="unif", adjustment = NULL, truncation=6000)
# fit a simple smooth of x and y
dsm.uni <- dsm(count~s(x, y), uni.mod, segdata, obsdata)
# second attempt using dummy_ddf
dsm.uni.dummy <- dsm(count~s(x, y), dummy_ddf(obsdata$object,width=6000), segdata, obsdata)
# Calculate the variance
# this will give a summary over the whole area in mexdolphins$preddata
mod.uni.var <- dsm_var_gam(dsm.uni, preddata, off.set=preddata$area)
summary(mod.uni.var)
mod.uni.var.dummy <- dsm_var_gam(dsm.uni.dummy, preddata, off.set=preddata$area)
summary(mod.uni.var.dummy)

Preferred behaviour

dsm_var_gam should recognise when a strip transect has been used and return simply the variability of the fitted gam because there is no uncertainty associated with the detection function model

[erex]

I edited print.summary.dsm.var to get around the troublesome lines of code (see below). I further pieced together some code to generate cell-specific CV values for subsequent plotting (when using a strip transect and dsm_var_gam see below:

### calling routine
load("E:/troubleshoot/Constanza/workspace.RData")
library(dsm)
dim(observation.data)[1]
dim(segment.data)[1]
dim(prediction.data)[1]
vis.gam(dsm.uni, plot.type = "contour", view = c("x","y"))

dsm.dummy <- dsm(count~s(x,y, k=15)+s(rug, k=10)+s(distcosta, k=10),family=tw(),
               dummy_ddf(observation.data$object,width=6000), 
               segment.data=segment.data, 
               observation.data=observation.data, method="REML")

dummy.var <- dsm_var_gam(dsm.dummy, prediction.data, 
                         off.set=prediction.data$area)
#
# before the final line, source the R function from the email:
source("E:/troubleshoot/Constanza/print.summary.dsm.var.ER.R")
# it replaces the function that causes the error when using the "dummy_ddf" object
#
print.summary.dsm.var.ER(summary(dummy.var))


# CV for each prediction grid cell
preddata.var <- split(prediction.data, 1:nrow(prediction.data))
var.forplot <- dsm_var_gam(dsm.dummy, pred.data=preddata.var,
                          off.set=prediction.data$area)
gridcell.cv <- sqrt(var.forplot$pred.var)/unlist(var.forplot$pred)
# I'm no good at plotting spatial data, but here is an object with x, y, CV for each cell
CV.by.cell <- cbind(prediction.data$x, prediction.data$y, gridcell.cv)
# Notice the size of the cell-by-cell CV values, very high uncertainty in the predictions
#   the smallest CV is 0.58, the median CV over all 409 cells is 2.37
summary(gridcell.cv)

modified print.summary.dsv.var

#' Print summary of density surface model variance object
#'
#' See [`summary.dsm.var`][summary.dsm.var] for information.
#'
#' @param x a summary of `dsm` variance object
#' @param \dots unspecified and unused arguments for S3 consistency
#' @return `NULL`
#' @export
#' @author David L. Miller
#' @seealso [`summary.dsm.var`][summary.dsm.var]
#' @keywords utility
print.summary.dsm.var.ER <- function(x, ...){
   if(x$bootstrap){
    cat("Summary of bootstrap uncertainty in a density surface model\n")
     if(!x$ds.uncertainty){
      cat("Detection function uncertainty incorporated using the delta method.\n\n")
     }else{
      cat("Detection function uncertainty incorporated into boostrap.\n\n")
    }
    cat("Boxplot coeff     :",x$boxplot.coef,"\n")
    cat("Replicates        :",x$n.boot,"\n")
    cat("Outliers          :",x$boot.outliers,"\n")
    cat("Infinites         :",x$boot.infinite,"\n")
    cat("NAs               :",x$boot.NA,"\n")
    cat("NaNs              :",x$boot.NaN,"\n")
    cat("Usable replicates : ",x$boot.usable,
        " (",100*(1-x$trim.prop),"%)\n",sep="")
      # Note that when we don't include detection function uncertainty in
    # our bootstrap then we need to do the log-Normal approx with the
    # CV calculated from the bootstrap, as the quantiles don't include
    # any of the uncertainty in the detection function.
    if(!x$ds.uncertainty){
      # asymptotic CI
      # this doesn't transform N, only se(N)
      # this probably should do the following:
      # lower =
      #  qlnorm(alpha/2, log(x$pred.est) - 0.5*log(x$cv^2+1),
      #         sqrt(log(x$cv^2+1)))
      # upper =
      #  qlnorm(1-alpha/2, log(x$pred.est) - 0.5*log(x$cv^2+1),
      #         sqrt(log(x$cv^2+1)))
      unconditional.cv.square <- x$cv^2
      asymp.ci.c.term <- exp(qnorm(1-x$alpha/2) *
                               sqrt(log(1+unconditional.cv.square)))
      asymp.tot <- c(x$pred.est / asymp.ci.c.term,
                     x$pred.est,
                     x$pred.est * asymp.ci.c.term)
      names(asymp.tot) <- c(paste0(x$alpha/2*100, "%"),
                            "Mean",
                            paste0((1-x$alpha/2)*100,"%"))
        cat("Approximate asymptotic bootstrap confidence interval:\n")
      print(asymp.tot)
      cat("(Using log-Normal approximation)\n")
      }else{
      # DSU percentile CI
      cat("\nPercentile bootstrap confidence interval and median:\n")
      print(x$quantiles)
    }
    }else if(!x$bootstrap){
    cat("Summary of uncertainty in a density surface model calculated\n")
    if(x$varprop){
      cat(" by variance propagation.\n")
      cat("\nProbability of detection in fitted model and variance model\n")
      print(x$saved$model.check)
    }else{
      cat(" analytically for GAM, with delta method\n")
      }
     cat("\n")
      ## calculate the CI around the abundance estimate
      # this doesn't transform N, only se(N)
    # this probably should do the following:
    # lower =
    #  qlnorm(alpha/2, log(x$pred.est) - 0.5*log(x$cv^2+1),
    #         sqrt(log(x$cv^2+1)))
    # upper =
    #  qlnorm(1-alpha/2, log(x$pred.est) - 0.5*log(x$cv^2+1),
    #         sqrt(log(x$cv^2+1)))
    unconditional.cv.square <- x$cv^2
    asymp.ci.c.term <- exp(qnorm(1-x$alpha/2) *
                             sqrt(log(1+unconditional.cv.square)))
    asymp.tot <- c(x$pred.est / asymp.ci.c.term,
                   x$pred.est,
                   x$pred.est * asymp.ci.c.term)
      names(asymp.tot) <- c(paste0(x$alpha/2*100, "%"),
                          "Mean",
                          paste0((1-x$alpha/2)*100,"%"))
      cat("Approximate asymptotic confidence interval:\n")
    print(asymp.tot)
    cat("(Using log-Normal approximation)\n")
    }
  if(class(x$saved$dsm.object$ddf)[1] != "fake_ddf") {
    cat("\n")
    cat("Point estimate                 :", x$pred.est,"\n")
    # print the individual CVs if we used the delta method
    if(x$bootstrap){
      if(!x$ds.uncertainty){
        if(!is.null(x$detfct.cv)){
          cat("CV of detection function       :",x$detfct.cv,"\n")
        }
        cat("CV from bootstrap              :", round(x$bootstrap.cv,4),"\n")
        cat("Total standard error           :", x$se,"\n")
        cat("Total coefficient of variation :", round(x$cv,4),"\n")
      }else{
        cat("Standard error                 :", x$se,"\n")
        cat("Coefficient of variation       :", round(x$cv,4),"\n")
      }
    }else{
      if(x$varprop | is.null(x$saved$dsm.object$ddf)){
        cat("Standard error                 :", x$se,"\n")
        cat("Coefficient of variation       :", round(x$cv,4),"\n")
      }else{
        if(!is.null(x$detfct.cv)){
          cat("CV of detection function       :",x$detfct.cv,"\n")
        }
        cat("CV from GAM                    :", round(x$gam.cv,4),"\n")
        cat("Total standard error           :", x$se,"\n")
        cat("Total coefficient of variation :", round(x$cv,4),"\n")
      }
    }
  }
  cat("Total standard error (only from GAM)          :", x$se, "\n")
  cat("Total CV (only from GAM)          :", round(x$se/x$pred.est,4), "\n")
  cat("\n")
    if(!is.null(x$subregions)){
    i <- 1
    for(subreg in x$subregions){
      cat("Subregion", i, "\n")
      print(subreg)
      cat("\n")
      i <- i+1
    }
  }
  invisible()
}