Why doesn't a binomial GLM successfully fix violated assumptions of simulated binomial data?

Question

I generated binomial data showing a logistic correlation to a predictor. I analysed such dataset with a generalised linear model which assumes a binomial distribution of residuals and a "logit" link function, but the residuals of such a model still fail to respect the model assumptions: the quantile-quantile plot shows that the residuals' distribution departs from the theoretical one, and residuals seem heteroscedastic. This happens regardless of sample size. Why is that?

Here is my example, simulated and analysed in R. To give the data some realism: let's imagine that I am interested in the probability of bird predation under different degrees of vegetation cover (expressed as percentage). For each level of vegetation cover I tracked 15 birds and recorded if they had been killed by a predator by the end of the study (data coded as 1 or 0).

# Generate data:
set.seed(666)
predictor <- c(0, 20, 40, 50, 70, 90)
y = 4 - 0.2 * predictor # linear combination of the variables
pr = 1/(1+exp(-y)) # pass y through an inv-logit function to get probability

# build dataset:

df <- data.frame(vegetation.cover = rep(predictor, 15), prob = rep(pr, 15))
df$predation.01 <- rbinom(length(df[,1]), 1, pr)
    
# df$predation.01 refers to the probability of a bird to be killed by a predator. 
# To model predation.01 vs vegetation.cover:

glm.01 <- glm(predation.01 ~ vegetation.cover,
    data= df,
    family=binomial(link="logit")
    )

These are the model diagnostic plots:

For comparison, these are the diagnostic plots for the same model fitted on a dataset with a replication of 1500 for each value of the predictor:

Here are the data together with the model predictions for n=15:

I added a pinch of artificial noise to the observations to avoid a perfect overlap of all the 0s and 1s in an attempt to aid the interpretation of the graph.

Answer 1

It may help you to read my answer to Interpretation of plot (glm.model). The plots that R is generating are not specific to a logistic regression. They are the plots for a standard OLS regression. For example, the qq-plot is not of the residuals vs a binomial distribution. It shows them relative to a normal distribution. However, there is no theoretical reason that the residuals from a logistic regression 'should' be normally distributed. On the other hand, the residuals from a logistic regression very much should be heteroscedastic, because the variance is a function of the mean. Your model seems fine, as far as I can tell.