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
|
require(robustbase)
set.seed(17)# reproducibility!
## to check:
## - for the empty model
summary(lmrob(Y ~ 0, coleman))
## - with and without an intercept in the model
summary(lmrob(Y ~ 1, coleman))
writeLines(sfm <- capture.output(
summary(lmrob(Y ~ ., coleman)))) # and this must be "identical":
sfm2 <- capture.output(summary(lmrob(Y ~ ., coleman, model=FALSE, x=FALSE, y=FALSE)))
iCall <- grep("lmrob.*coleman", sfm)# the only line that differs
stopifnot(sfm[-iCall] == sfm2[-iCall])
## w/o intercept:
summary(lmrob(Y ~ . - 1, coleman, model=FALSE, x=FALSE, y=FALSE))
## - when prior-weights are included
wts <- c(rep(0.05, 10), rep(2, 10))
summary(lmrob(Y ~ . - 1, coleman, model=FALSE, x=FALSE, y=FALSE,
weights = wts))
## - should work for object with NA in the coefficients, and
## - should work for object with NA in the observations --> both in ./NAcoef.R
## check equality with lm() for classical model
test <- function(formula, data,
items=c("coefficients", "residuals", "df", "scale",
"r.squared", "adj.r.squared", "weights"),
tol = 1e-4, ...)
{
lmrCtrl <- lmrob.control(psi = "hampel", tuning.psi = c(1000, 2000, 3000),
method="SMDM", ...)
sc <- summary(lm (formula, data))
sr <- summary(lmrob(formula, data, control= lmrCtrl))
names(sc)[names(sc) == "sigma"] <- "scale"
if(sc$df[1] == 0 && getRversion() <= "3.5.1" && as.numeric(R.version$`svn rev`) < 74993)
## in the past, lm() returned logical empty matrix
storage.mode(sc$coefficients) <- "double"
ret <- all.equal(sc[items], sr[items], tolerance=tol)
if (!isTRUE(ret)) {
print(sr)
for (i in seq_along(items)) {
print(sc[items[i]])
print(sr[items[i]])
}
print(ret)
stop(sprintf("all.equal(sc[items], sr[items], tol.. = %g) are not all TRUE",
tol))
}
ret
}
set.seed(101)
test(Y ~ 0, coleman, c("residuals", "df", "coefficients",
"r.squared", "adj.r.squared"), tol=1e-10)
test(Y ~ 1, coleman, tol = 2e-4)
test(Y ~ ., coleman, tol = 4e-4)
test(Y ~ . - 1, coleman, tol = 4e-4)
|
