Fatal error when using dsm_var_prop

[erex]

Error arises on line 11 of varprop_check

dsm_var_prop produces a reasonable standard error for the Gomex dolphins, but the last line of dsm_var_prop is to return a result including a call to varprop_check that generates the error.

Attempt to manufacture MRE from the Gomex data. Requires artificially creating a Beaufort value at the segment level (requirement of using dsm_var_prop for variance estimation).

library(dsm)
library(ggplot2)
load("mexdolphins-extra.rda") # must be downloaded from distanceexamples website
data(mexdolphins)
library(rgdal)
library(maptools)
library(plyr)
# tell R that the survey.area object is currently in lat/long
proj4string(survey.area) <- CRS("+proj=longlat +datum=WGS84")
# proj 4 string
# using http://spatialreference.org/ref/esri/north-america-lambert-conformal-conic/
lcc_proj4 <- CRS("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")
# project using LCC
survey.area <- spTransform(survey.area, CRSobj=lcc_proj4)
# simplify the object
survey.area <- data.frame(survey.area@polygons[[1]]@Polygons[[1]]@coords)
names(survey.area) <- c("x", "y")
# given the argument fill (the covariate vector to use as the fill) and a name,
# return a geom_polygon object
# fill must be in the same order as the polygon data
grid_plot_obj <- function(fill, name, sp){
# what was the data supplied?
names(fill) <- NULL
row.names(fill) <- NULL
data <- data.frame(fill)
names(data) <- name
spdf <- SpatialPolygonsDataFrame(sp, data)
spdf@data$id <- rownames(spdf@data)
spdf.points <- fortify(spdf, region="id")
spdf.df <- join(spdf.points, spdf@data, by="id")
# seems to store the x/y even when projected as labelled as
# "long" and "lat"
spdf.df$x <- spdf.df$long
spdf.df$y <- spdf.df$lat
geom_polygon(aes_string(x="x",y="y",fill=name, group="group"), data=spdf.df)
}

library(Distance)
detfc.hr.null <- ds(distdata, max(distdata$distance), key="hr", adjustment=NULL)
gof_ds(detfc.hr.null)

detfc.hr.beau<-ds(distdata, max(distdata$distance), formula=~beaufort,
                      key="hr", adjustment=NULL)
gof_ds(detfc.hr.beau)

segdata$beaufort <- rpois(nrow(segdata), 1.2)
segdata$beaufort <- ifelse(segdata$beaufort<5, segdata$beaufort, 4)
dsm.beau.xy <- dsm(count~s(x,y), detfc.hr.beau, segdata, obsdata, method="REML")
summary(dsm.beau.xy)
dsm.xy.pred <- predict(dsm.beau.xy, preddata, preddata$area)

preddata.var <- split(preddata, 1:nrow(preddata))
dsm.xy.var <- dsm_var_gam(dsm.beau.xy, pred.data=preddata.var,
                          off.set=preddata$area)
summary(dsm.xy.var)  # var_gam is happy

dsm.xy.varprop <- dsm_var_prop(dsm.beau.xy, pred.data=preddata.var,
                          off.set=preddata$area)  #var_prop is not happy
summary(dsm.xy.varprop)

I don't know the purpose of varprop_check so I don't know if simply eliminating the call to the function in the last line of dsm_var_prop is a viable solution to the problem.

Other inconsistencies: object names in the list returned by dsm_var_prop do not correspond with object names used in summary.dsm_varprop; e.g. first element of list is pred.var, but summary thinks there is an element var

[VLucet]

Hello, @crentmeister, @StewartResearch and I are working to apply your methods to our dataset (thank you for the very nice documentation and package by the way), and we have been running into that same issue. We are interested in using the variance propagation method so we don't have to assume independence between the distance function estimation and the dsm, so this fix is important for us.

I have made an attempt to fix that bug. I haven't touched the core of the code, I just tried to make the objects conform to the expected format from my limited understanding of the code, making things into lists or not depending. It seems to work (using the same MRE above) but it feels like a hack.

Below is the output I get:

Summary of uncertainty in a density surface model calculated
 by variance propagation.

Probability of detection in fitted model and variance model
[[1]]
  beaufort Fitted.model Fitted.model.se Refitted.model
1      1.2    0.7037521       0.2687439      0.6611248
2      3.0    0.5840859       0.2427054      0.5877802
3      4.8    0.4687253       0.3683304      0.5160213


Approximate asymptotic confidence interval:
      2.5%       Mean      97.5% 
  470.5135  3742.5847 29769.4775 
(Using log-Normal approximation)

Point estimate                 : 3742.585 
Standard error                 : 5375.688 
Coefficient of variation       : 1.4364 

It doesn't seem to match the output in your Supplementary Materials, but I think the model is different (simpler?) in that MRE (but I cant be sure the code is the same as the SM seem to use deprecated versions of those functions).

FWIW, the tests run fine on the PR branch:

==> devtools::test()

ℹ Testing dsm
This is dsm 2.3.3.9001
Built: NA
✔ | F W S  OK | Context
✔ |         2 | Do GLMs work [0.4s]                                                           
✔ |        18 | test inputs [4.0s]                                                            
✔ |         6 | Mexico pantropical dolphin data [10.7s]                                       
✔ |         5 | prediction [0.5s]                                                             
✔ |         3 | GAM variance [1.8s]                                                           

══ Results ═══════════════════════════════════════════════════════════════════════════════════
Duration: 17.5 s

[ FAIL 0 | WARN 0 | SKIP 0 | PASS 34 ]

[lenthomas]

Thanks @VLucet, @crentmeister and @StewartResearch - very helpful. I took an initial look and concluded that it's going to take quite a bit of digging to get up to speed on the objects that are passed in to this function -- I'd like to understand how this issue arose. I'm away for a week but will have a look when I get back. In the meantime, I assume you folks are fine using your forked version, and are making further progress with your work.