diags.lm: Standard diagnostics for lm objects

In gmonette/spida: Collection of miscellaneous functions for mixed models etc. prepared for SPIDA 2009+ (development version)

Description

Standard diagnostics for lm objects

Usage

diags.lm(x, y, ..., ask, labels = names(residuals(x)), showlabs = text)

Arguments

x
y
...
ask
labels
showlabs

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, ..., ask, labels = names(residuals(x)), showlabs = text) 
{
    if (!missing(ask)) {
        op <- par(ask = ask)
        on.exit(par(op))
    }
    form <- formula(x)
    f <- predict(x)
    r <- residuals(x)
    nams <- names(r)
    if (!missing(labels)) {
        nams <- names(residuals(x))
        if (length(nams) != length(labels)) 
            labels <- labels[nams]
    }
    ret <- NULL
    if (missing(y)) {
        y <- f + r
        yname <- deparse(form[[2]])
    }
    else yname <- deparse(substitute(y))
    fname <- paste("Fitted:", deparse(form[[3]]), collapse = " ")
    plot(f, y, xlab = fname, ylab = yname, main = "Dependent var. vs. Predicted", 
        ...)
    abline(0, 1, lty = 1)
    lines(supsmu(f, y))
    showlabs(f, y, labels, ...)
    lmi <- lm.influence(x)
    hat <- lmi$hat
    sigma <- lmi$sigma
    mm <- scale(model.matrix(x), scale = F)
    mp <- predict(x, type = "terms")
    comp.res <- mp + r
    plot(f, abs(r), xlab = fname, ylab = deparse(substitute(abs(resid(x)))), 
        main = "Absolute Residual vs. Predicted", ...)
    showlabs(f, abs(r), labels, ...)
    zq <- qqnorm(r, main = "Normal Quantile Plot", ylab = "Residual", 
        sub = fname)
    qqline(r)
    showlabs(zq, labels, ...)
    n <- length(r)
    r.o <- sort(r)
    half <- (n + 1)/2
    if (n%%2 == 1) {
        med <- r.o[half]
        below <- med - r.o[half:1]
        above <- r.o[half:n] - med
    }
    else {
        med <- sum(r.o[c(half, half + 1)])/2
        below <- med - r.o[(n/2):1]
        above <- r.o[(n/2 + 1):n] - med
    }
    opt <- par(pty = "s")
    ran <- range(c(below, above))
    plot(below, above, main = "Symmetry plot of residuals", xlab = "Distance below median", 
        ylab = "Distance above median", xlim = ran, ylim = ran)
    abline(0, 1, lty = 2)
    par(opt)
    std.r <- r/(sigma * sqrt(1 - hat))
    plot(hat, std.r, xlab = "Leverage (hat)", ylab = yname, sub = fname, 
        main = "Studentized residual vs. Leverage", ...)
    showlabs(hat, std.r, labels, ...)
    nams <- dimnames(lmi$coefficients)[[1]]
    pairs(lmi$coefficients)
    pairs(lmi$coefficients, panel = function(x, y, nams) {
        points(x, y)
        text(x, y, nams)
    }, nams = nams)
    invisible(0)
  }

gmonette/spida documentation built on May 17, 2019, 7:25 a.m.