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.


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

# 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.

Main residual output

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

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.

Residual against Model Prediction Residual against fitted(m1)

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?

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.