1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
## MM: more following of plot.lm() : ~/R/D/r-devel/R/src/library/stats/R/plot.lm.R
plot.lmrob <-
function (x, which = 1:5,
caption = c("Standardized residuals vs. Robust Distances",
"Normal Q-Q vs. Residuals", "Response vs. Fitted Values",
"Residuals vs. Fitted Values" ,
"Sqrt of abs(Residuals) vs. Fitted Values"),
panel = if(add.smooth) panel.smooth else points,
sub.caption = deparse(x$call), main = "",
compute.MD = TRUE, # maybe (n < 1000 && p < 20)
ask = prod(par("mfcol")) < length(which) && dev.interactive(),
id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
label.pos = c(4,2), qqline = TRUE, add.smooth = getOption("add.smooth"),
..., p = 0.025)
{
if (!inherits(x, "lmrob"))
stop("Use only with 'lmrob' objects")
if (!is.numeric(which) || any(which < 1) || any(which > 5))
stop("'which' must be in 1:5")
show <- rep(FALSE, 5)
show[which] <- TRUE
r <- residuals(x)
n <- length(r)
sr <- r/x$scale
yh <- fitted(x)
if (is.null(id.n))
id.n <- 0
else {
id.n <- as.integer(id.n)
if(id.n < 0L || id.n > n)
stop(gettextf("'id.n' must be in {1,..,%d}", n), domain = NA)
}
if(id.n > 0L) { ## label the largest residuals
if(is.null(labels.id))
labels.id <- paste(1L:n)
iid <- 1L:id.n
show.r <- sort.list(abs(r), decreasing = TRUE)[iid]
## if(any(show[2L:3L]))
## show.rs <- sort.list(abs(rs), decreasing = TRUE)[iid]
text.id <- function(x, y, ind, adj.x = TRUE) {
labpos <-
if(adj.x) label.pos[1+as.numeric(x > mean(range(x)))] else 3
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]) {
if(is.null(x[['MD']]) && compute.MD) {
message("recomputing robust Mahalanobis distances")
x$MD <- ## need to recompute
robMD(x = if(!is.null(x[['x']])) x$x else
if(!is.null(x[['model']])) model.matrix(x, x$model)
else stop("need 'model' or 'x' component for robust Mahalanobis distances"),
intercept = attr(x$terms,"intercept"),
wqr = x$qr)
## try to "cache" them with the object
.ge <- .GlobalEnv
if(identical(parent.frame(), .ge) &&
exists((cnx <- as.character(match.call()[["x"]])), .ge)) {
assign(cnx, x, envir = .ge)
message("saving the robust distances 'MD' as part of ", sQuote(cnx))
}
}
if(!is.null(xD <- x[['MD']])) {
if (p < 0 || p > 1)
stop ("Tolerance range must be between 0% to 100%")
else chi <- sqrt( qchisq(p = 1-p, df = x$rank) )
ylim <- range(sr, na.rm=TRUE)
if(id.n > 0) ylim <- extendrange(r = ylim, f = 0.08)
plot(xD, xlab = "Robust Distances",
sr, ylab = "Robust Standardized residuals", ylim=ylim,
main = main, type = "n", ...)
panel(xD, sr, ...)
mtext(caption[1], 3, 0.25)
if (one.fig)
title(sub = sub.caption, ...)
if(id.n > 0) {
y.id <- sr[show.r]
y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
text.id(xD[show.r], y.id, show.r)
}
abline(h = c(2.5,-2.5), lty = 3)
abline(v = chi, lty = 3)
}
}
if (show[2L]) { ## Normal
qq <- qqnorm(r, ylab = "Residuals", main = main,...)
if(qqline) qqline(r, lty = 3, col = "gray50")
mtext(caption[2], 3, 0.25)
if (one.fig)
title(sub = sub.caption, ...)
if(id.n > 0)
text.id(qq$x[show.r], qq$y[show.r], show.r)
}
if (show[3]) {
y <- if(!is.null(x[['model']])) model.response(x$model) else yh + r
m1 <- min(yh,y)
m2 <- max(yh,y)
plot(yh, y, xlab = "Fitted Values", ylab = "Response",
xlim = c(m1,m2), ylim = c(m1,m2), main = main, type = "n", ...)
panel(yh, y, ...)
mtext(caption[3], 3, 0.25)
if (one.fig)
title(sub = sub.caption, ...)
if(id.n > 0)
text.id(yh[show.r], y[show.r], show.r)
abline(a = 0,b = 1)
}
if (show[4]) {
plot(yh, r, xlab = "Fitted Values", ylab = "Residuals",
main = main, type = "n", ...)
panel(yh, r, ...)
mtext(caption[4], 3, 0.25)
if (one.fig)
title(sub = sub.caption, ...)
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)
}
abline(h = c(2.5*x$scale,0,-2.5*x$scale), lty = 3)
}
if (show[5]) {
sqrtabsr <- sqrt(abs(r))
plot(yh, sqrtabsr, xlab = "Fitted Values", ylab = "Sqrt of abs(Residuals)",
main = main, type = "n", ...)
panel(yh, sqrtabsr, ...)
mtext(caption[5], 3, 0.25)
if (one.fig)
title(sub = sub.caption, ...)
if(id.n > 0)
text.id(yh[show.r], sqrtabsr[show.r], show.r)
}
invisible()
}
|
