brms compatibility#523
Conversation
…il in spatialAutocorrelation plot so that margins are returned to previous settings
- adding getPearsonResiduals.brmsfit and improving testDispersion to detect absence of Pearson residuals in brms models and return NA for tests with refit=T AND type="PearsonChisq" - adding info in help file of getData to alert users that it will search data frames in the environment.
this functions doesn't interfere with others, but to keep it standardized I created the compatibility for brms because brms::fixef() returns a matrix (mean, se, q2.5, q97.5) and the standard for other packages is to return a numeric vector. When refit=T, the simulateResiduals was creating a weird matrix in simulationOutput$refittedFixedEffects. Now it is ok.
using the function getPearsonResiduals for the refitted Pearsons to 1) use the function we created for that to avoid errors/warnings from residuals from model's package. This way we guarantee NAs in brms models, for example.
Reason: if I have getPearsonResiduals.brmsfit and don't export it - R CMD CHECK will complain with a note on S3 methods. If I export it, it will appear in the help of getPearsonResiduals, but it doesn't seem elegant/nice to have a function there that just throws an error. So, I decided to implement the error for brmsfit models in getPearsonResiduals.default.
…hartig/DHARMa into v0.5.0_brmsCompatibility
…l 01 / bernoulli and k/n data
…r, testData is removed and created anew for phylolm); adjusted vignette for bayesians for brms compatibility
"brms and Stan models are unique because they compile C++ code during the test phase. Even if BH is a dependency of your dependency, the compiler sometimes fails to link it unless it's explicitly present in the environment's library."
adding rstan to suggests instead of Enhances
adding BH and StanHeaders to suggestions to see if the brms testthat works
adding RcppEigen in extra-packages
deleting the packages I added to `suggests`, to test if it will still work in GH actions
melina-leite
force-pushed
the
v0.5.0_brmsCompatibility
branch
from
acaa3a5 to
4745ac2
Compare
About the GH Action fails and the solutionWith brms, in order to make the testthat run, we had to add the bold code in R-CMD-check.yaml: uses: r-lib/actions/setup-r-dependencies@v2 However, there was a suggestion (from AI) to also add the same packages to the 'Suggests' section of the DESCRIPTION FILE. In the end, I kept the latest commit without this suggestion. I'm not sure the consequences of this to CRAN submission. |
melina-leite
force-pushed
the
v0.5.0_brmsCompatibility
branch
from
802e7c5 to
b8f35e1
Compare
| message("Note that the Chi2 test on Pearson residuals is biased for MIXED models towards underdispersion. Tests with alternative = two.sided or less are therefore not reliable. If you have random effects in your model, we recommend to test only with alternative = 'greater', i.e. test for overdispersion, or else use the DHARMa default tests which are unbiased. See help for details.")} | ||
|
|
||
|
|
||
| #rp <- getPearsonResiduals(model) |
There was a problem hiding this comment.
Hi @cosminawerneke - I wonder if we really need this error message here or if it would not be sufficient to throw the error in getPearsonResiduals, and then people can see the traceback?
In the way you coded it it is of course slightly more informative / cleaner, but also higher maintenance.
There was a problem hiding this comment.
Hi @florianhartig,
I added this here (for type="Pearson") as it was also in the (type = "DHARMa" & refit=T). Cause the error wasn't clear when having no Pearson residuals in this case.
I'm not sure about moving it to getPearsonResiduals, because, as I see it, we would have to repeat it in every getPearsonResiduals.xx function (there aren't that many, but still repeating code).
Let me know if you have a more elegant way of doing it! :)
There was a problem hiding this comment.
hmm ... but it still feels a bit weird - getPearsonResiduals is our own function, so I feel we should be able to call it and rely on it providing sensible output without having to catch errors.
Or, in other words, getPearsonResiduals should be implemented for all sensible model classes, and specific errors should be thrown. For the default function, we could catch a general error.
I would still change this. If this code is used somewhere else in the package, we should change it as well, it may been a remnant from the time where I didn't have a getPearsonResiduals function
There was a problem hiding this comment.
@florianhartig, changes made:
-
I moved the error message to getPearsonResiduals.default and getPearsonResiduals.gam
-
I had to suppress the error message in simulateResiduals() with refit=T to avoid repeating the error message nSim times. Then, I throw just ONE warning message at the end of simulateResiduals(refit=T).
Here is a working example with mgcv:
devtools::load_all()
data <- createData()
# DOESN'T HAVE PEARSON - shouldn't work
library(mgcv)
mod2 <- gam(observedResponse ~ Environment1,
data=data, family = ziP())
summary(mod2)
mgcv:::residuals.gam(mod2, type="scaled.pearson") # shouldn't work
# error message of mgcv is weird, it searches for res, but because the family doesn't have a formula for Pearson, it simply says res wasn't found. (could be more informative, but it is not)
getPearsonResiduals(mod2) # ok, not working: DHARMa message + mgcv error message
testDispersion(mod2, type="Pearson") # ok the same
resis <- simulateResiduals(mod2)
getPearsonResiduals(resis$fittedModel) # ok the same
testDispersion(resis, type="DHARMa") # dharma test works ok
resis2 <- simulateResiduals(mod2, refit=T) # only ONE error message in the end
testDispersion(resis2, type="DHARMa") # ok the same
There was a problem hiding this comment.
Hi @florianhartig, I implemented the changes we discussed today (commit 70be55f. Please have a look at it.
There was a problem hiding this comment.
Moreover, I had to change a bit the testthat for some packages that don't have the Pearson (brms, GLMMadaptive, phylolm) to avoid the error message to stop the R CMD check.
|
Hi @cosminawerneke, looks fine for me except for the few questions I have made directly in the code. One general question: I didn't see anything that tests if the brms model that is used is suitable for DHARMa - what will happen if a person uses e.g. a multinomial model etc. ? |
|
Hi @florianhartig, as we agreed on in our last meeting about brms compatibility, DHARMa does not test if a brms model can reliably be checked or not. We just give a warning every time a brms model is detected (see compatibility.R line 45), and then it's up to the user to decide if their model is suitable to be used with DHARMa |
|
@melina-leite, I would still clean up this error message, otherwise ready to be merged. |
- moving Pearson residuals error messagem from testDisperstion to getPearsonResiduals() - suppressing error message of not getting pearson in simulateResiduals when refit=T - some gitignore from new version rstudio "positai"
- adding the getPearsonResiduals.brms(), although it throws an error automatically because Pearson residuals are deprecated for brms. I added a note about it in the help file of the wrapper. - deleted the message about not having Pearson residuals when doing refit=T in simulateResiduals(). - Cleaning unnecessary code in testDispersion.
given that the getPearsonResiduals.brms function will always give an error, I created a runEverything function excluding the Pearson residuals dispersion tests in testModelTypes. Used this for brms, phylolm/phyloglm, GLMMadaptive. Otherwise the tests wouldn't work
3 participants
Compatibility for simple brms models (i.e. models which could also be fit with glmmTMB; no multinomial, multi-response or structural equation models)