jranke / chemCal
1
calplot <- function(object, 
2
  xlim = c("auto", "auto"), ylim = c("auto", "auto"), 
3
  xlab = "Concentration", ylab = "Response",
4
  legend_x = "auto", 
5
  alpha = 0.05, varfunc = NULL)
6
{
7 0
  UseMethod("calplot")
8
}
9

10
calplot.default <- function(object, 
11
  xlim = c("auto","auto"), ylim = c("auto","auto"), 
12
  xlab = "Concentration", ylab = "Response",
13
  legend_x = "auto",
14
  alpha=0.05, varfunc = NULL)
15
{
16 0
  stop("Calibration plots only implemented for univariate lm objects.")
17
}
18

19
calplot.lm <- function(object,
20
  xlim = c("auto","auto"), ylim = c("auto","auto"), 
21
  xlab = "Concentration", ylab = "Response",
22
  legend_x = "auto", 
23
  alpha=0.05, varfunc = NULL)
24
{
25 0
  if (length(object$coef) > 2)
26 0
    stop("More than one independent variable in your model - not implemented")
27

28 0
  if (alpha <= 0 | alpha >= 1)
29 0
    stop("Alpha should be between 0 and 1 (exclusive)")
30

31 0
  m <- object
32 0
  level <- 1 - alpha
33 0
  y <- m$model[[1]]
34 0
  x <- m$model[[2]]
35 0
  if (xlim[1] == "auto") xlim[1] <- 0
36 0
  if (xlim[2] == "auto") xlim[2] <- max(x)
37 0
  xlim <- as.numeric(xlim)
38 0
  newdata <- list(
39 0
    x = seq(from = xlim[[1]], to = xlim[[2]], length=250))
40 0
  names(newdata) <- names(m$model)[[2]]
41 0
  if (is.null(varfunc)) {
42 0
    varfunc <- if (length(m$weights)) {
43 0
        function(variable) mean(m$weights)
44 0
      } else function(variable) rep(1,250)
45
  }
46 0
  pred.lim <- predict(m, newdata, interval = "prediction",
47 0
    level=level, weights.newdata = varfunc(m))
48 0
  conf.lim <- predict(m, newdata, interval = "confidence",
49 0
    level=level)
50 0
  yrange.auto <- range(c(0,pred.lim))
51 0
  if (ylim[1] == "auto") ylim[1] <- yrange.auto[1]
52 0
  if (ylim[2] == "auto") ylim[2] <- yrange.auto[2]
53 0
  if (legend_x[1] == "auto") legend_x <- min(object$model[[2]])
54 0
  plot(1,
55 0
    type = "n", 
56 0
    xlab = xlab,
57 0
    ylab = ylab,
58 0
    xlim = as.numeric(xlim),
59 0
    ylim = as.numeric(ylim)
60
  )
61 0
  points(x,y, pch = 21, bg = "yellow")
62 0
  matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), 
63 0
    col = c("black", "red", "red"))
64 0
  if (length(object$weights) > 0) {
65 0
    legend(min(x), 
66 0
      max(pred.lim, na.rm = TRUE), 
67 0
      legend = c("Fitted Line", "Confidence Bands"),
68 0
      lty = c(1, 3), 
69 0
      lwd = 2, 
70 0
      col = c("black", "green4"), 
71 0
      horiz = FALSE, cex = 0.9, bg = "gray95")
72
  } else {
73 0
  matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), 
74 0
    col = c("black", "green4", "green4"))
75 0
  legend(legend_x,
76 0
    max(pred.lim, na.rm = TRUE), 
77 0
    legend = c("Fitted Line", "Confidence Bands", 
78 0
        "Prediction Bands"), 
79 0
    lty = c(1, 3, 4), 
80 0
    lwd = 2, 
81 0
    col = c("black", "green4", "red"), 
82 0
    horiz = FALSE, cex = 0.9, bg = "gray95")
83
  }
84
}

Read our documentation on viewing source code .

Loading