"plot.lm" <-
function (x, which = 1:4, caption = c("Residuals vs Fitted", 
                              "Normal Q-Q", "Scale-Location",
                              "Cook's distance", "Residuals vs Leverage"), 
            panel = points, sub.caption = deparse(x$call), main = "", 
            ask = prod(par("mfcol")) < length(which) && dev.interactive(), 
            ..., id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
            cook.levels=c(0.5, 1.0), label.pos=c(4,2))
{
  if (!inherits(x, "lm")) 
    stop("Use only with 'lm' objects")
  if (!is.numeric(which) || any(which < 1) || any(which > 5)) 
    stop("`which' must be in 1:5")
  isGlm <- inherits(x, "glm")
  show <- rep(FALSE, 5)
  show[which] <- TRUE
  r <- residuals(x)
  yh <- predict(x)
  w <- weights(x)
  if (!is.null(w)) {
    wind <- w != 0
    r <- r[wind]
    yh <- yh[wind]
    w <- w[wind]
    labels.id <- labels.id[wind]
  }
  n <- length(r)
  if (any(show[2:5])) {
    s <- if (inherits(x, "rlm")) 
      x$s
    else sqrt(deviance(x)/df.residual(x))
    hii <- lm.influence(x, do.coef = FALSE)$hat
    if (any(show[4:5])) {
      cook <- if (isGlm)cooks.distance(x)
      else cooks.distance(x, sd = s, res = r)
    }
  }
  if (any(show[c(2:3,5)])) {
    ylab23 <- if (isGlm) 
      "Std. deviance resid."
    else "Standardized residuals"
    r.w <- if (is.null(w)) 
      r
    else sqrt(w) * r
    rs <- r.w/(s * sqrt(1 - hii))
  }
  if (any(show[c(1, 3)])) 
    l.fit <- if (isGlm) 
      "Predicted values"
    else "Fitted values"
  if (is.null(id.n)) 
    id.n <- 0
  else {
    id.n <- as.integer(id.n)
    if (id.n < 0 || id.n > n) 
      stop("`id.n' must be in {1,..,", n, "}")
  }
  if (id.n > 0) {
    if (is.null(labels.id)) 
      labels.id <- paste(1:n)
    iid <- 1:id.n
    show.r <- sort.list(abs(r), decreasing = TRUE)[iid]
    if (any(show[2:3])) 
      show.rs <- sort.list(abs(rs), decreasing = TRUE)[iid]
    text.id <- function(x, y, ind, adj.x = FALSE){
      midx <- mean(range(x))
      labpos <- if (!adj.x) 3
      else
        label.pos[1+as.numeric(x>midx)]
      text(x, y, labels.id[ind], cex = cex.id, xpd = TRUE, pos=labpos,
           offset=0.25)
    }
  }
  one.fig <- prod(par("mfcol")) == 1
  if (ask) {
    op <- par(ask = TRUE)
    on.exit(par(op))
  }
  if (show[1]) {
    ylim <- range(r, na.rm = TRUE)
    if (id.n > 0) 
      ylim <- ylim + c(-1, 1) * 0.08 * diff(ylim)
    plot(yh, r, xlab = l.fit, ylab = "Residuals", main = main, 
         ylim = ylim, type = "n", ...)
    panel(yh, r, ...)
    if (one.fig) 
      title(sub = sub.caption, ...)
    mtext(caption[1], 3, 0.25)
    if (id.n > 0) {
      y.id <- r[show.r]
      y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
      text.id(yh[show.r], y.id, show.r, adj.x = TRUE)
    }
    abline(h = 0, lty = 3, col = "gray")
  }
  if (show[2]) {
    ylim <- range(rs, na.rm = TRUE)
    ylim[2] <- ylim[2] + diff(ylim) * 0.075
    qq <- qqnorm(rs, main = main, ylab = ylab23, ylim = ylim, 
                 ...)
    if (one.fig) 
      title(sub = sub.caption, ...)
    mtext(caption[2], 3, 0.25)
    if (id.n > 0) 
      text.id(qq$x[show.rs], qq$y[show.rs], show.rs, adj.x = TRUE)
  }
  if (show[3]) {
    sqrtabsr <- sqrt(abs(rs))
    ylim <- c(0, max(sqrtabsr, na.rm = TRUE))
    yl <- as.expression(substitute(sqrt(abs(YL)), list(YL = as.name(ylab23))))
    yhn0 <- if (is.null(w)) 
      yh
    else yh[w != 0]
    plot(yhn0, sqrtabsr, xlab = l.fit, ylab = yl, main = main, 
         ylim = ylim, type = "n", ...)
    panel(yhn0, sqrtabsr, ...)
    if (one.fig) 
      title(sub = sub.caption, ...)
    mtext(caption[3], 3, 0.25)
    if (id.n > 0) 
      text.id(yhn0[show.rs], sqrtabsr[show.rs], show.rs, 
              adj.x = TRUE)
  }
  if (show[4]) {
    if (id.n > 0) {
      show.r <- order(-cook)[iid]
      ymx <- cook[show.r[1]] * 1.075
    }
    else ymx <- max(cook)
    plot(cook, type = "h", ylim = c(0, ymx), main = main, 
         xlab = "Obs. number", ylab = "Cook's distance", ...)
    if (one.fig) 
      title(sub = sub.caption, ...)
    mtext(caption[4], 3, 0.25)
    if (id.n > 0) 
      text.id(show.r, cook[show.r], show.r, adj.x=FALSE)
  }
  if (show[5]) {    
    ylim <- range(rs, na.rm = TRUE)
    hatval <- hatvalues(x)
    if (id.n > 0) {
      ylim <- ylim + c(-1, 1) * 0.08 * diff(ylim)
      show.r <- order(-cook)[iid]
    }
    plot(hatval, rs, ylim = ylim, main = main, 
         xlab = "Leverages", ylab = ylab23,
         type="n", ...)
    panel(hatval, rs, ...)
    if (one.fig) 
      title(sub = sub.caption, ...)
    p <- length(coef(x))
    for(crit in cook.levels){
      curve(sqrt(crit*p*(1-x)/x), lty=2, add=T)
      curve(-sqrt(crit*p*(1-x)/x), lty=2, add=T)
    }
    xmax <- par()$usr[2]
    ymult <- sqrt(p*(1-xmax)/xmax)
    aty <- c(-sqrt(rev(cook.levels))*ymult, sqrt(cook.levels)*ymult)
    axis(4, at=aty, labels=paste(c(rev(cook.levels), cook.levels)),
         mgp=c(.25,.25,0), las=2, tck=0, cex.axis=cex.id)
    mtext(caption[5], 3, 0.25)
    if (id.n > 0) {
      y.id <- rs[show.r]
      y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
      text(hatval[show.r], y.id, paste(show.r), pos=2, cex=cex.id, offset=0.25)
  }
  }
  if (!one.fig && par("oma")[3] >= 1) 
    mtext(sub.caption, outer = TRUE, cex = 1.25)
  invisible()
}