How can I cure non-normality of residuals in logistic regression? I am running a logistic regression against the "Default" dataset of ISLR.
I am using the performance package in R (available on CRAN) to test the goodness of my model.
MWE:
rm(list=ls())
indf <- read.csv("https://pastebin.com/raw/suXpddVR")
indf$default <- as.factor(indf$default)
indf$student <- as.factor(indf$student)
mod <- glm(default ~ student+balance+income,
           data=indf,
           family="binomial")
print(summary(mod))

library(performance)
performance::check_model(mod)

Diagnostic plot is as follows:

As you can see from the plot in the upper right the residuals don't follow a normal distribution.
How can I cure that?
I have tried putting the squares of the variables in the model, but they are non-significant and the problem remains.
 A: I haven't checked what the performance package calculates exactly, but these are likely deviance or Pearson residuals. Both behave pretty similar and are based on the idea to scale the residual deviation to the expected variance.
Most GLM distributions (such as Poisson or binomial) will be asymptotically normal with some dependency of the variance on the mean, and in this case, deviance / Pearson residuals will become normal. So, for large counts (Poisson) or a large number of trials (binomial), it makes sense to check those residuals for normality.
Unfortunately, those requirements are rarely met in practice.  In your case, I expect that you have 0/1 data, and in which case this won't work at all.
A possible solution would be to aggregate data points. A far better solution, however, is to calculate randomised quantile residuals, which have proper distributional properties across the entire spectrum of possible outcomes. I recommend using the DHARMa package (disclaimer: I am the developer). Read the vignette, which explains how this works, and note that there is a separate section on 0/1 logistic residuals.
