Broken function: glmResidPlot

Claude.ai, after some prompting, suggests:

glm_diagnostics <- function(glmobj, plots = c(1, 2, 3, 4), layout = NULL) {
    if (!inherits(glmobj, "glm")) {
        stop("Object has not been obtained from fitting a glm")
    }

    famobj <- glmobj$family$family
    mu.trans <- switch(famobj,
        binomial = 2 * asin(sqrt(glmobj$fitted.values)),
        gaussian = glmobj$fitted.values,
        Gamma = 2 * log(glmobj$fitted.values),
        poisson = 2 * sqrt(glmobj$fitted.values),
        stop("Residuals for glms using ", famobj, " are not available.")
    )

    xlabel <- switch(famobj,
        binomial = expression(paste("2 arcsine(", sqrt(widehat(mu)), ")")),
        gaussian = expression(widehat(mu)),
        Gamma = expression(paste("2 log", widehat(mu))),
        poisson = expression(paste("2 ", "sqrt(widehat(mu))"))
    )

    resi <- stats::residuals(glmobj, type = "deviance")
    hatvals <- stats::influence(glmobj)$hat
    idNum <- seq_along(resi)
    hatOne <- (hatvals >= 1)
    resi <- resi[!hatOne]
    hatvals <- hatvals[!hatOne]
    mu.trans <- mu.trans[!hatOne]
    idNum <- idNum[!hatOne]
    resi <- resi/sqrt(1 - hatvals)

    plot_functions <- list(
        function() {
            plot(mu.trans, resi, xlab = xlabel, ylab = "Standardised deviance residual", 
                pch = 16, bty = "l", las = 1, main = "Residuals vs. transformed fitted means")
            graphics::abline(h = 0, col = "gray50")
            low <- stats::lowess(mu.trans, resi, iter = 0)
            graphics::lines(low, col = "red", lwd = 2)
        },
        function() {
            plot(idNum, resi, xlab = "Index number", ylab = "Standardised residual", 
                pch = 16, bty = "l", las = 1, main = "Index plot of residuals")
            graphics::abline(h = 0, col = "gray50")
        },
        function() {
            plot(idNum, resi^2, col = 1 + (resi < 0), pch = 16, xlab = "Index number", 
                ylab = "Squared standardised residual", bty = "l", 
                las = 1, main = "Index of squared residuals")
        },
        function() {
            stats::qqnorm(resi, main = "Normal Q-Q plot", las = 1, bty = "l", pch = 16)
            stats::qqline(resi)
        }
    )

    plots <- plots[plots %in% 1:4]
    n_plots <- length(plots)

    if (is.null(layout)) {
        layout <- switch(as.character(n_plots),
                         "1" = c(1, 1),
                         "2" = c(1, 2),
                         "3" = c(2, 2),
                         "4" = c(2, 2))
    }

    graphics::par(mfrow = layout)

    for (i in plots) {
        plot_functions[[i]]()
    }
}

Test:

options(device = function(...){
    png(...)
    dev.control("enable")
 }, webr.plot.new = FALSE)

library(M348)

y <- c(10,8,9,2,5,0)
x <- c(1,2,3,4,5,7)
glmfit <- glm(y~x, family="poisson")

And then:

all plots: glm_diagnostics(glmfit)
select plots: glm_diagnostics(glmfit, plots = c(1, 3))

ouseful-demos / jupyterlite-m348-demo

Broken function: glmResidPlot #3