Add support for piecewiseSEM in DHARMa

[florianhartig]

I have been thinking about this, see comments in https://github.com/jslefche/piecewiseSEM/issues/197

Perspective for this request: Currently, the package does not allow to simulate from the fitted model, which prevents an easy integration into DHARMa. Support could be added if the package developers implement a simulate function, see here

Interim solution for users that, for some reason, have to urgently use DHARMa with this package: take your fitted model, create a simulate function for this model structure yourself, and then use createDHARMa (see help), this will allow most options of the package to be run. See further comments on support of new packages in here. Note that the vignette has some further comments / examples on creating custom simulation functions and reading them into DHARMa

[florianhartig]

Here a piece of code that calculates DHARMa residuals for a fitted psem object

library(piecewiseSEM)

mod = psem(
  lm(rich ~ distance + elev + abiotic + age + hetero + firesev + cover, data =  keeley),
  lm(firesev ~ elev + age + cover, data =  keeley), 
  lm(cover ~ age + elev + abiotic, data =  keeley)
)

summary(mod)
plot(mod)


library(DHARMa)
residuals2 <- function(fittedPSEM, type = "DHARMa"){
  if(type == "standard"){
    return(residuals(fittedPSEM))
  }
  if(type == "DHARMa"){
    sim = list()
    res = list()
    for(i in 1:(length(fittedPSEM) - 1)){
      sim[[i]] <- simulateResiduals(fittedPSEM[[i]])
      res[[i]] = sim[[i]]$scaledResiduals
    }
    res = matrix(unlist(res),ncol=length(res),byrow=F)
    return(list(residuals = res, DHARMaSimulations = sim))
  }
}

res = residuals2(mod)
res$residuals
# plots and tests for model 1
plot(res$DHARMaSimulations[[1]])
testResiduals(res$DHARMaSimulations[[1]])
# in principle, it would be possible to plot residuals also against a particular path, or test for homogeneity along this path