Open florianhartig opened 4 years ago
We just discussed this in some length last week in a working group at ISEC 2020. In principle, if you can simulate from your fitted model, you can create residuals for each entry of the community matrix with the createDHARMa
function (just note that you will have to ensure that the simulate function returns a vector, DHARMa cannot handle matrix responses), and based on this, you can still use all DHARMa plots, so I would do
aggregateResiduals
to group residuals by site and species, and check that as well with all standard plots / checksThere are probably a few more things that one could do to check the specific structure of a jSDM, and that was the content of our discussion. If something happens on this front, I will update this issue.
I tried the univariate approach (i.e. dropping the dimensions and treating the JSDM outputs univariate) but it doesn't produce right residual plots. I then tried a multivariate ecdf and it seems to produce correct scaled residuals (perhaps this is different for LVM-JSDMs where we assume a univariate likelihood):
library(sjSDM) # devel branch
library(DHARMa)
library(mltools)
library(data.table)
## simulate community:
com = simulate_SDM(env = 2L, species = 7L, sites = 100L)
## fit model:
model = sjSDM(Y = com$response,
env = linear(com$env_weights, ~.),
iter = 50L,
family = binomial("probit"))
# increase iter for your own data
sims = simulate(model, nsim = 200L)
dim(sims)
sp = 1
dim(sims)
sim =
createDHARMa(t(array(sims, c(200, 100*7))),
c(com$response),
c(predict(model)),
integerResponse = TRUE)
plot(sim)
vs =
sapply(1:100, function(i) {
empirical_cdf(as.data.table(sims[,i,]) + matrix(runif(200*7, -0.5, 0.5), 200, 7),
ubounds = as.data.table(com$response[i,]+runif(7,-0.5,0.5)))$CDF
})
sim$scaledResiduals = c(vs)
sim$observedResponse = c(com$response)
sim$fittedPredictedResponse = c(predict(model))
sim$nObs = 700
plot(sim)
So the solution for JSDM support could be to implement the multivariate ecdf in DHARMa.
Hi Max,
thanks - 2 questions:
a) in your multivariate code, when you use matrix(runif(200*7, -0.5, 0.5), 200, 7), this is the old idea of DHARMa randomisation, in all new versions I'm using the PIT randomisation per default ... is there a specific reason why you switched to the old randomisation procedure?
b) On a first glance, I don't really get what you mean by "treating the outputs multivariate" - you to the ecdf per site (sapply), but then in sapply, you create one residual per species, right? I think I have to run it myself to see what's in the values, but just to get the idea ...
Add-on note: in principle we could aggregate residuals per site (over species) or per species (over sites), or show all residuals, as in the internal structure ...
Question from a user