DHARMa integration

[florianhartig]

Hi Jon,

I'm the developer of DHARMa and incidentally just attending a workshop on SEMs, which made me think about the possibility to integrate DHARMa with piecewiseSEM.

I have looked at the structure of your package, and I see two routes one could go

1. Integrating DHARMa into piecewiseSEM: most of the models that are supported by piecewiseSEM could also be checked with DHARMa (except nlme). So, it would be possible to add switch in piecewiseSEM and return DHARMa residuals for each model, or p-values for DHARMa omnibus residual tests that could be superimposed on the DAG or displayed in piecewiseSEM.

2. Integrating piecewiseSEM into DHARMa: to support piecewiseSEM, minimum requirements are described here https://github.com/florianhartig/DHARMa/wiki/Adding-new-R-packages-to-DHARMA

Option 2 would allow to create omnibus for a tested SEM from DHARMa, e.g. global dispersion of spatial autocorrelation tests. Given the nature of the package, however, it seems to me more logical to consider each link individually (because averaging over them might mask problems), which would suggest to me that it would be more logical to go via route 1. I think what could added easily is implement

a) the ability to return DHARMa instead of normal residuals per link b) the option to impose the results of testResiduals (3 general tests) or testUniformity (tests just for distribution) on the links in the dag, e.g. in the summary() or in the plot() function, to alert users to connections that are problematic

[florianhartig]

Hi, I just played a bit around with that - here an example of how one could modify the residuals function in piecewiseSEM to support DHARMa ... maybe one should check in the loop if the fitted models are supported by DHARMa, and fall back to standard otherwise.

Add-on: obviously, the model doesn't make a lot of sense, and also DHARMa doesn't add much for lms, but just as a proof of principle.

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

[jslefche]

This is clever and timely. I would be interested in integrating more robust tests of assumptions in the package –looping over the list seems the most applicable since there are no formal omnibus tests across multiple models. Its good you opened an issue as that will keep it on my radar

---

Jonathan S. Lefcheck, Ph.D. Tennenbaum Coordinating Scientist MarineGEO: https://marinegeo.si.edu/ Smithsonian Institution Phone: +1 (443) 482-2443 www.jonlefcheck.net

From: Florian Hartigmailto:[email protected] Sent: Tuesday, December 17, 2019 7:20 AM To: jslefche/piecewiseSEMmailto:[email protected] Cc: Subscribedmailto:[email protected] Subject: Re: [jslefche/piecewiseSEM] DHARMa integration (#197)

External Email - Exercise Caution

Hi, I just played a bit around with that - here an example of how one could modify the residuals function in piecewiseSEM to support DHARMa ... maybe one should check in the loop if the fitted models are supported by DHARMa, and fall back to standard otherwise.

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

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fjslefche%2FpiecewiseSEM%2Fissues%2F197%3Femail_source%3Dnotifications%26email_token%3DAAR4AVZ2N25RVLL62LRLUJLQZC7YDA5CNFSM4J3ZGVPKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHCGGQQ%23issuecomment-566518594&data=02%7C01%7Clefcheckj%40si.edu%7C57692a5a903f4ce0428808d782eb7ad7%7C989b5e2a14e44efe93b78cdd5fc5d11c%7C0%7C0%7C637121820201441666&sdata=OvadqVztUJQblFdJ%2Bfl%2Fk21DNmeR%2FD7dx2IoiWxTZgM%3D&reserved=0, or unsubscribehttps://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAR4AV4VQZLPH7KFMN334GTQZC7YDANCNFSM4J3ZGVPA&data=02%7C01%7Clefcheckj%40si.edu%7C57692a5a903f4ce0428808d782eb7ad7%7C989b5e2a14e44efe93b78cdd5fc5d11c%7C0%7C0%7C637121820201451665&sdata=ljGVU1AA7Udpf3zFO7sjhYWSyCVJiu8ISaM6k%2FUnqgE%3D&reserved=0.

[jebyrnes]

Yup - I've just lapplied over objects in the past. This approach should work great.

[florianhartig]

I guess you would have to decide if you just want to return the scaled residuals, or a list with residuals and the DHARMa objects (as I did here).

There forme makes it compatible to the normal residuals, but one looses the possibility to run the tests implemented in DHARMa on the returned residuals.