Skip to content
Prev 275694 / 398506 Next

show.ci='FALSE' ignored in simple.lm

Hi:

A trip to package sos reveals that simple.lm() is in the UsingR
package. Looking at the code of this function, it plots the (x, y)
pairs and the fitted least squares line without an option to suppress
the plot. Here's a slight hack of the function; it adds a new argument
plot with default action TRUE; after the linear model is fit, an if
statement tests to see if the plot is desired; if FALSE, it simply
returns the (print method of the) fitted object.

simple.lm2 <- function (x, y, plot = TRUE, show.residuals = FALSE,
          show.ci = FALSE, conf.level = 0.95, pred = FALSE)
{
    op <- par()
    ord <- order(x)
    x <- x[ord]
    y <- y[ord]
    tmp.lm <- lm(y ~ x)
    if(!plot) return(tmp.lm) else {

    if (show.residuals) {
        par(mfrow = c(2, 2))
    }
    plot(x, y)
    abline(tmp.lm)

    if (show.ci) {
        xvals <- seq(min(x), max(y), length = 15)
        curve(predict(tmp.lm, data.frame(x = x), level = conf.level,
            interval = "confidence")[, 3], add = TRUE)
        curve(predict(tmp.lm, data.frame(x = x), level = conf.level,
            interval = "confidence")[, 2], add = TRUE)
        curve(predict(tmp.lm, data.frame(x = x), level = conf.level,
            interval = "prediction")[, 3], lty = 3, add = TRUE)
        curve(predict(tmp.lm, data.frame(x = x), level = conf.level,
            interval = "prediction")[, 2], lty = 3, add = TRUE)
    }
    coeffs <- floor(tmp.lm$coeff * 100 + 0.5)/100
    plusorminus <- c("+")
    if (coeffs[1] < 0)
        plusorminus <- c("")
    title(paste("y = ", coeffs[2], "x ", plusorminus, coeffs[1]))
    if (show.residuals) {
        Fitted <- fitted.values(tmp.lm)
        Residuals <- residuals(tmp.lm)
        plot(Fitted, Residuals)
        abline(h = 0)
        title("Residuals vs. fitted")
        hist(Residuals, main = "hist of residuals")
        qqnorm(Residuals, main = "normal plot of residuals")
        qqline(Residuals)
    }
    if (pred) {
        print(predict(tmp.lm, data.frame(x = pred)))
    }
      }
    tmp.lm
}


#####
# Simple test:

plot(x <- 1:10)
y <- 5*x + rnorm(10,0,1)

tmp1 <- simple.lm2(x, y, plot = FALSE)
summary(tmp1)
# Produces the plot
simple.lm2(x, y)



HTH,
Dennis
On Thu, Oct 27, 2011 at 7:35 AM, David Winsemius <dwinsemius at comcast.net> wrote: