I would like to find out what variables influence/explain an efficiency score for an invasive species control method. As the score is defined on (0,1] and as some observations were not independent, I wanted to use one-inflated beta regression mixed-models using
I thus ran the following code to test this function with a simple model:
library(gamlss) mod1 <- gamlss(formula = efficiency~distance+re(random = ~1|manager_id), nu.formula = ~distance+re(random = ~1|manager_id), data = eff, family = BEOI())
I used the same formula for the mu and nu parts of the model as I had no reason to suspect that the probabilities of having an outcome at (0,1) or at 1 should be driven by different variables.
At this point, I tried to inspect the residuals and assess the goodness-of-fit of my model and that's when things started to get ugly. First, the residuals had a rather strange look:
Second, the Chi-square goodness-of-fit test I performed (following this tutorial) indicated a lack of fit and yet, the computed pseudo-R2 was incredibly high...:
dat.resid <- sum(resid(mod1, type = "weighted", what = "mu")^2) 1 - pchisq(dat.resid, mod1$df.resid)  0.000628797 Rsq(mod1)  0.8490156
... clearly suggesting overfitting (
gamlss::Rsq()computes the generalised R-squared of Nagelkerke (1991) which is, unless I'm mistaken, not adapted to mixed-models but still, this is weird):
My questions are thus:
- What could be the cause of such apparent model misspecification? With n = 85, my sample size should be large enough to use this type of model, shouldn't it? Is the
gamlsspackage known for being unreliable?
- Could someone guide me as to how to correctly perform diagnostic checks on one-inflated beta regression mixed-models? According to Mark White (see here and here), beta regression models are relatively new and thus lack consensual diagnostic methods. However, his answers are now 3 years old and perhaps someone developed since then some robust procedure to assess the quality of this type of model (possibly using R)?
- What could be the cause of the discrepancy between the Pearson’s Chi-square test result and the too-good-to-be-true fit displayed by the R2 and the observed vs. fitted plot?
As you can probably see, I'm not very experienced in statistics so any correction, documentation and constructive criticism is welcome. Thanks in advance.
If someone wants to reproduce my results, here are my data:
library(RCurl) library(dplyr) x <- RCurl::getURL("https://raw.githubusercontent.com/mrelnoob/jk.dusz.tarping/master/mydata/erad.csv") read.csv(text = x) %>% mutate(manager_id = as.factor(manager_id)) %>% # There are other factors in the dataset, but it doesn't matter here mutate(efficiency = efficiency/10) -> eff # For my Y to be [0,1]