Residual diagnostics in DHARMa for multilevel logistic regression

Question

I'm running a multilevel logistic regression, and have been trying to look at residual diagnostics using the DHARMa package. My data is copy-pastable from here and the following code should run everything.

library(lme4)
library(DHARMa)

m1 <- glmer("error ~ 1 + year + categorisation + statistic + (1 | journalID)",
                 data = data,
                 family = binomial(link = "logit"))

simulationOutput <- simulateResiduals(fittedModel = m1)

# Main plot function from DHARMa, which gives 
# Left: a qq-plot to detect overall deviations from the expected distribution
# Right: a plot of the residuals against the rank-transformed model predictions
plot(simulationOutput)

# Plotting standardized residuals against predictors
plotResiduals(simulationOutput, data$year, xlab = "Year", main=NULL)
plotResiduals(simulationOutput, data$statistic, xlab = "Statistic Type", main=NULL)
plotResiduals(simulationOutput, data$categorisation, xlab = "Category 1 or Category 2", main=NULL)
plotResiduals(simulationOutput, data$journalID, xlab = "Journal ID", main=NULL)

# Plotting standardized residuals against the predicted value
plotResiduals(simulationOutput, main=NULL)
plotResiduals(simulationOutput, fitted(m1), xlab = "fitted(m1)")

The following is the main residual plot function.

There follow plots of standardized residuals against my continuous predictor Year, a nominal factor with 4 levels (Statistic), and against the binary predictor Category.

There follows the default DHARMa plot of residuals against the predicted values. I also tried to plot the residuals against fitted(m1), excepting that they would be the same thing. However, they aren't.

My questions are:

1. Why can't I plot the residuals against the predictor I modelled as a random effect (journalID). When I try to do it I get the error message Error in plot.window(...) : need finite 'xlim' values.
2. Why is it that my two different ways of standardized residuals against the predicted value, plotResiduals(simulationOutput, main=NULL), and plotResiduals(simulationOutput, fitted(m1), xlab = "fitted(m1)"), give slightly different results. Aren't they the same thing?
3. The plot for my predictor statistic lists the levels of that factor as 1-4, while the actual names of the levels of that factor are "F", "chi","t","r". How can I work out which number maps onto which level of the factor?
4. Do my residuals look OK? Is there anything further I should be testing in relation to them?

Answer 1

I'm the developer of DHARMa. Below answers to your questions, but note for the future that it would be preferable to ask questions about syntax problems on the DHARMa development site here.

1. Because in your data, data$journalID is coded as a character, if you code as.factor it works

2. No they are not, as there are different ways to make predictions (conditional / unconditional) for a mixed model. The way it's coded in DHARMa is highly preferable, as plotting against the default of fitted() can produce spurious patterns. See explanation here.

3. Use as.numeric(factor). It's true that I should maybe change this to displaying the factor names though.

4. From the plots, nothing springs to mind, but this is not uncommon for the binomial. Often useful to aggregate residuals for the 0/1 binomial. See vignette, there is also a section on the binomial model.

p.s.: I have now included an automatic fix for 1) and 3). This will be included in DHARMa 0.3.4. Until the, you could get it if you install the development version from here.