florianhartig / DHARMa

Diagnostics for HierArchical Regession Models
http://florianhartig.github.io/DHARMa/
200 stars 21 forks source link

Error when using simulateResiduals() on a phyr object (class communityPGLMM) #377

Closed emyoung90 closed 1 year ago

emyoung90 commented 1 year ago

Hi here,

I'm having issues getting DHARMa to play nicely with phyr with respect to the simulateResiduals() function. Per Daijiang Li's vignette, the following code should produce the residuals I need for standard diagnostic plots:

library(phyr)
library(ape)
library(dplyr)
library(DHARMa)
data("oldfield") #contained within the phyr package

mod_bayes <- phyr::pglmm(pres ~ disturbance + (1 | sp__) + (1 | site) + 
                         (disturbance | sp__) + (1 | sp__@site), 
                         data = oldfield$data, 
                         cov_ranef = list(sp = oldfield$phy),
                         family = "binomial",
                         bayes = TRUE,
                         prior = "pc.prior.auto")

resids <- DHARMa::simulateResiduals(mod_bayes, plot = FALSE)

However, despite using the same code and data, I'm repeatedly getting the following error: "Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'". I've made sure all packages etc are up-to-date, my environment is clean, etc.

I'm assuming that is because the model output, which is a list of 42 objects, contains some items that are "NULL". However, attempting to simply remove the null values (which I suspect is a bad idea anyway) gives me other errors. Predictably, I'm getting the same errors when I try this with my own data, so it's not just a quirk of the oldfield dataset.

Overall, I'm stumped. Is this a bug? Or are DHARMa and phyr no longer compatible?

florianhartig commented 1 year ago

Hi El,

your example doesn't run for me either. It suspect the problem is the Bayesian fit, when I run an example from me:

data("oldfield")

data = oldfield$data %>%
  dplyr::filter(disturbance == 1)

mod_disturbed <- phyr::pglmm(pres ~ (1 | sp__) , 
                             data = oldfield$data %>%
                               dplyr::filter(disturbance == 1), 
                             cov_ranef = list(sp = oldfield$phy),
                             family = "binomial")

res <- simulateResiduals(mod_disturbed)
plot(res) # works

this works without problems. The reason for the failure is that there seems to be no simulate() function implemented for models with with INLA, so this is a problem that needs to be solved on the phyr / INLA side.

I'll link this to the general phyr discussion and will close this issue!