Incorrect McKelvey & Zavoina pseudo-R2

[imbhe]

I believe the implementation of the McKelvey & Zavoina pseudo-R2 is incorrect. The numbers seem to be far too low. The variance calculated as np.var(y_pred) should actually be the variance of the logit-probabilities. Compare, e.g., to the implementation in R using DescTools::PseudoR2:
y.hat <- predict(x, type = "link").
s2 <- switch(link, probit = 1, logit = pi^2/3, NA)
sse <- sum((y.hat - mean(y.hat))^2)
res["McKelveyZavoina"] <- sse/(n * s2 + sse)