Two bugs in predict.dsm() for "density" response

[dfifield]

There are two problems with predict.dsm() when response is "density":

See reprex below for examples.

1. predict(dsm_model) returns 0

   When newdata and off.set are both NULL in the call to predict.dsm(), off.set never gets set properly before it is used in the code chunk at lines 78-87.

   Not that you need suggestions for how to fix it, but I feel bad just reporting bugs and not offering up any solutions, so

   Changing line 27 to:

       newdata <- object$data
       if(is.null(off.set)){
            off.set <- 1
       }
   

   fixes this problem.

2. predict(dsm_model, type = "lpmatrix") returns an error:

   Error in if (type == "response" & se.fit) {: argument is of length zero

   I stumbled across this while trying to get vis.concurvity() to work. This happens at line 80 in predict.dsm.R because se.fit is being tested when NULL. The simplest solution appears to be to change line 80 to:

   }else if(type=="response" & (!is.null(se.fit) && se.fit)){

   so that it's safe when type == "lpmatrix"

### Reprex follows
library(Distance)
#> Loading required package: mrds
#> This is mrds 2.2.0
#> Built: R 3.6.1; ; 2019-09-28 03:24:13 UTC; windows
#> 
#> Attaching package: 'Distance'
#> The following object is masked from 'package:mrds':
#> 
#>     create.bins
library(dsm)
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
#> Loading required package: numDeriv
#> This is dsm 2.2.18
#> Built: R 3.6.1; ; 2019-10-16 15:41:41 UTC; windows

sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> Platform: i386-w64-mingw32/i386 (32-bit)
#> Running under: Windows 10 x64 (build 17134)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dsm_2.2.18          numDeriv_2016.8-1.1 mgcv_1.8-28        
#> [4] nlme_3.1-140        Distance_0.9.8      mrds_2.2.0         
#> 
#> loaded via a namespace (and not attached):
#>  [1] minqa_1.2.4        BB_2014.10-1       ggplot2_3.2.1     
#>  [4] tidyselect_0.2.5   lattice_0.20-38    dfoptim_2018.2-1  
#>  [7] Rcgmin_2013-2.21   pkgconfig_2.0.1    compiler_3.6.1    
#> [10] stringr_1.4.0      parallel_3.6.1     Rcpp_1.0.2        
#> [13] optimx_2013.8.7    plyr_1.8.4         highr_0.6         
#> [16] tools_3.6.1        truncnorm_1.0-8    rmarkdown_1.16    
#> [19] Rvmmin_2018-4.17   R6_2.2.2           purrr_0.3.2       
#> [22] knitr_1.25         scales_0.5.0       svUnit_0.7-12     
#> [25] assertthat_0.2.0   digest_0.6.15      evaluate_0.14     
#> [28] gtable_0.2.0       setRNG_2013.9-1    Matrix_1.2-14     
#> [31] stringi_1.1.7      optextras_2016-8.8 htmltools_0.4.0   
#> [34] munsell_0.4.3      grid_3.6.1         glue_1.3.1        
#> [37] colorspace_1.3-2   rlang_0.4.0        quadprog_1.5-5    
#> [40] Rsolnp_1.16        magrittr_1.5       splines_3.6.1     
#> [43] lazyeval_0.2.1     yaml_2.2.0         crayon_1.3.4      
#> [46] xfun_0.9           dplyr_0.8.3        pillar_1.4.2      
#> [49] statmod_1.4.32     tibble_2.1.3       ucminf_1.1-4
data(mexdolphins)
hr.model <- ds(distdata, max(distdata$distance), key = "hr", adjustment = NULL)
#> Fitting hazard-rate key function
#> Key only model: not constraining for monotonicity.
#> AIC= 841.253
#> No survey area information supplied, only estimating detection function.

# fit a simple smooth of x and y to densities
mod1 <- dsm(density~s(x,y)+s(depth), hr.model, segdata, obsdata)
predict(mod1)
#> numeric(0)
predict(mod1, type = "lpmatrix")
#> Error in if (type == "response" & se.fit) {: argument is of length zero

[erex]

Thanks Dave. Dave Miller is out of the office this week. He'll be the fellow who'll follow up on making the repairs you suggest.

[dfifield]

Thanks Eric. In the meantime I've implemented the suggested changes in my fork. If it's helpful, I can issue a pull request for Dave. However, I expect it'll take him longer to read this issue report than it will to just fix it himself directly. :-)

[dill]

Not clear if you get the notifications for the commits @dfifield, but I've pushed a fix to this, do let me know if this works for you!

Thanks again for the thorough bug report (as usual)!

[dfifield]

@dill I had a quick look but I'm bit confused.

For example, in the following code, all of the predict calls return exactly the same thing.

library(Distance)
library(dsm)

data(mexdolphins)
hr.model <- ds(distdata, max(distdata$distance), key = "hr", adjustment = NULL)

# fit a simple smooth of x and y to densities
mod1 <- dsm(density~s(x,y)+s(depth), hr.model, segdata, obsdata)
predict(mod1, newdata = segdata)
predict(mod1, newdata = segdata, off.set = 1)
predict(mod1, newdata = segdata, off.set = 2)
predict(mod1, newdata = segdata, off.set = 20)
segdata$off.set <-  5
predict(mod1, newdata = segdata)

It looks like newdata$off.set is being nicely computed in your code but I can't see where it's being used.

Is this intentional? Am I missing something?

[dill]

Sorry Dave, that code is painful and I messed-up. The docos previously said that the values in newdata should be used over those in off.set. I've changed this to be the other way around, since that makes more sense. Give this one a whirl and see if that's better.

I think really I need to just remove the off.set argument, and make people have a column in the prediction data with the area in it. Something for the next version...

[dfifield]

So sorry to keep bugging you on this one Dave!

I'm still having two problems with the latest version.

Problem 1
With a density response model, I get the same output from predict(model, newdata, off.set = x) regardless of what I set x to. This doesn't happen with abundance response models. Note that it doesn't matter if off.set is supplied in newdata or as an arg to predict().

Is this the intended behavior? I think it isn't, but I'm a bit confused at this point so I want to make sure I'm not off base. Should predict.dsm() be returning the same thing regardless of what off.set is for density response models?

Suggested solution
Recently (two commits ago) you removed the code at the bottom of predict.dsm() that formerly multiplied the return value from predict.gam() by the offset (for density models only). Does that multiplication need to be re-instated?

Here's a reprex.

library(Distance)
#> Loading required package: mrds
#> This is mrds 2.2.0
#> Built: R 3.6.1; ; 2019-09-28 03:24:13 UTC; windows
#> 
#> Attaching package: 'Distance'
#> The following object is masked from 'package:mrds':
#> 
#>     create.bins
library(dsm)
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
#> Loading required package: numDeriv
#> This is dsm 2.2.18
#> Built: R 3.6.1; ; 2019-10-22 15:30:43 UTC; windows

data(mexdolphins)
hr.model <- ds(distdata, max(distdata$distance), key = "hr", adjustment = NULL)
#> Fitting hazard-rate key function
#> Key only model: not constraining for monotonicity.
#> AIC= 841.253
#> No survey area information supplied, only estimating detection function.

#  First an abundance model which works fine with predict.
mod1 <- dsm(abundance.est~s(x,y)+s(depth), hr.model, segdata, obsdata)

# just so I don't have to create a new dataframe I'll re-use segdata as newdata
newdata <- segdata

# Let's compare output of predict() with two different offsets
x <- predict(mod1, newdata, off.set = 2)
y <- predict(mod1, newdata, off.set = 4)

# Predictions in y should be twice those in x. Are they?
head(y/x) # Equals 2. Correct.
#> 1 2 3 4 5 6 
#> 2 2 2 2 2 2

# Now fit a density response model. Predict with this model does not give what I expect.
mod2 <- dsm(density~s(x,y)+s(depth), hr.model, segdata, obsdata)

x <- predict(mod2, newdata, off.set = 2)
y <- predict(mod2, newdata, off.set = 4)

# Should y = 2x just like in the abundance model above?
head(y/x) # Equals 1. Incorrect??
#> 1 2 3 4 5 6 
#> 1 1 1 1 1 1


# It doesn't matter if off.set is supplied in newdata or as an arg to predict()
newdata$off.set <- 2
x <- predict(mod2, newdata)

newdata$off.set <- 4
y <- predict(mod2, newdata)

# Same result.
head(y/x) # Equals 1. Incorrect??
#> 1 2 3 4 5 6 
#> 1 1 1 1 1 1

Problem Two
As of the most recent commit (e9af1ae), if newdata is not supplied, then the off.set argument to predict.dsm() is being ignored (actually over-written, see below).

Suggested solution
I think there's an else clause missing for the if statement on lines 32-34 of predict.dsm().

Here's the chunk of code in question. Note that if off.set is supplied then newdata$off.set is set to the value of off.set, but is then immediately overwritten.

# if there is no newdata, then use the data stored in the model object
  if(is.null(newdata)){
    newdata <- object$data

    if(!is.null(off.set)){
      newdata$off.set <- off.set
    }
    newdata$off.set <- object$family$linkfun(newdata$segment.area)
  }else{

I have to admit, I don't quite understand the bit with the ...$linkfun and whether it should be applied to off.set or not, but that's another question....

I'm guessing you intended it to be something like:

# if there is no newdata, then use the data stored in the model object
  if(is.null(newdata)){
    newdata <- object$data

    if(!is.null(off.set)){
      newdata$off.set <- off.set
    } else {
      newdata$off.set <- object$family$linkfun(newdata$segment.area)
    }
  }else{

Here's a reprex to demonstrate the problem:

library(Distance)
#> Loading required package: mrds
#> This is mrds 2.2.0
#> Built: R 3.6.1; ; 2019-09-28 03:24:13 UTC; windows
#> 
#> Attaching package: 'Distance'
#> The following object is masked from 'package:mrds':
#> 
#>     create.bins
library(dsm)
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
#> Loading required package: numDeriv
#> This is dsm 2.2.18
#> Built: R 3.6.1; ; 2019-10-22 15:30:43 UTC; windows

data(mexdolphins)
hr.model <- ds(distdata, max(distdata$distance), key = "hr", adjustment = NULL)
#> Fitting hazard-rate key function
#> Key only model: not constraining for monotonicity.
#> AIC= 841.253
#> No survey area information supplied, only estimating detection function.

#  First an abundance model which works fine with predict.
mod1 <- dsm(abundance.est~s(x,y)+s(depth), hr.model, segdata, obsdata)

# Note these all return the same thing.
head(predict(mod1, off.set = 1))
#>           1           2           3           4           5           6 
#> 0.005324244 0.006400757 0.006854569 0.007908613 0.009628911 0.010993812
head(predict(mod1, off.set = 2))
#>           1           2           3           4           5           6 
#> 0.005324244 0.006400757 0.006854569 0.007908613 0.009628911 0.010993812
head(predict(mod1, off.set = 3))
#>           1           2           3           4           5           6 
#> 0.005324244 0.006400757 0.006854569 0.007908613 0.009628911 0.010993812

Thanks so much for all your efforts on this and the package in general. It's greatly appreciated!

[dfifield]

Thanks for the update Dave!

There's one small issue left that breaks vis.concurvity() for density models only. This was my original "Problem 2" above.

vis.concurvity() eventually causes predict.dsm(object, type = "lpmatrix") to be executed, which for density models, gives the following error in predict.dsm():

Error in if (type == "response" & se.fit) {: argument is of length zero

This happens at line 79 in predict.dsm.R because type == "lpmatrix" and se.fit is NULL.

Changing this line to:

}else if(type=="response" & (!is.null(se.fit) && se.fit)){

fixed the problem for me.

Here's a reprex to demonstrate the problem:

library(Distance)
#> Loading required package: mrds
#> This is mrds 2.2.0
#> Built: R 3.6.1; ; 2019-09-28 03:24:13 UTC; windows
#> 
#> Attaching package: 'Distance'
#> The following object is masked from 'package:mrds':
#> 
#>     create.bins
library(dsm)
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
#> Loading required package: numDeriv
#> This is dsm 2.2.18
#> Built: R 3.6.1; ; 2019-10-24 12:57:54 UTC; windows

data(mexdolphins)
hr.model <- ds(distdata, max(distdata$distance), key = "hr", adjustment = NULL)
#> Fitting hazard-rate key function
#> Key only model: not constraining for monotonicity.
#> AIC= 841.253
#> No survey area information supplied, only estimating detection function.

#  First an abundance model which works fine with predict.
mod1 <- dsm(density~s(x,y)+s(depth), hr.model, segdata, obsdata)

# Error
vis.concurvity(mod1)
#> Error in if (type == "response" & se.fit) {: argument is of length zero
predict(mod1, type = "lpmatrix")
#> Error in if (type == "response" & se.fit) {: argument is of length zero

[dill]

Apologies Dave, I was half-way through the previous message when I got pulled into a series of meetings. I've just pushed some changes that should fix what you mention above. Thanks for your patience!

[dfifield]

Oops, I jumped the gun there - so thanks for YOUR patience! Next time I'll sit on my hands and be patient until I'm told you're done :-)

Thanks for all your fixes - it's working for the test cases I've tried.

Closing #37...