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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
|
R Under development (unstable) (2026-02-01 r89366) -- "Unsuffered Consequences"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ## Test handling of NA coefficients / singular fits
> ## also check:
> ## -- what would have to be done if class "lm" was added.
> ## -- general compatibility to class lm.
>
> require(robustbase)
Loading required package: robustbase
> options(digits = 5, width = 111) -> op # -> higher chance of platform independence
>
> ## generate simple example data (almost as in ./weights.R )
> data <- expand.grid(x1=letters[1:3], x2=LETTERS[1:3], rep=1:3)
> set.seed(1)
> data$y <- rnorm(nrow(data))
> ## drop all combinations of one interaction:
> data <- subset(data, x1 != 'c' | (x2 != 'B' & x2 != 'C'))
> ## add collinear variables
> data$x3 <- rnorm(nrow(data))
> data$x4 <- rnorm(nrow(data))
> data$x5 <- data$x3 + data$x4
> ## add some NA terms
> data$y[1] <- NA
> data$x4[2:3] <- NA ## to test anova
>
> ## Classical models start with 'cm', robust just with 'rm' (or just 'm'):
> cm0 <- lm (y ~ x1*x2 + x3, data)
> cm1 <- lm (y ~ x1*x2 + x3 + x4 + x5, data)
> set.seed(2)
> rm1 <- lmrob(y ~ x1*x2 + x3 + x4 + x5, data)
> m3 <- lmrob(y ~ x1*x2 + x3 + x4, data) # same column space as rm1
> rm0 <- lmrob(y ~ x1*x2 + x3, data)
>
> ## clean version of rm1 (to check predict)
> data2 <- data.frame(y=data$y[-(1:3)], rm1$x[,!is.na(rm1$coef)])
> set.seed(2)
> rm1c <- lmrob(y ~ x1b + x1c + x2B + x2C + x3 + x4 + x1b:x2B + x1b:x2C, data2)
>
> ## add class lm to rm1 (for now)
> class(rm1) <- c(class(rm1), "lm")
> class(rm0) <- c(class(rm0), "lm")
>
> ## the full matrix (data) should be returned by model matrix (frame)
> stopifnot(all.equal(model.matrix(cm1), model.matrix(rm1)),
+ all.equal(model.frame (cm1), model.frame (rm1)))
> ## qr decomposition should be for the full data and pivots identical lm result
> qr.cm1 <- qr(cm1)$qr
> qr.rm1 <- rm1$qr$qr
> stopifnot(NCOL(qr.rm1) == NCOL(qr.cm1),
+ NROW(qr.rm1) == NROW(qr.cm1),
+ length(rm1$qr$qraux) == length(qr(cm1)$qraux),
+ all.equal(rm1$qr$pivot, qr(cm1)$pivot),
+ all.equal(dimnames(qr.rm1),dimnames(qr.cm1)))
> ## the alias function should return the same result
> stopifnot(all.equal(alias(cm1), alias(rm1)))
>
> ####
> ## these helper functions should print NAs for the dropped coefficients
> print(rm1)
Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data)
\--> method = "MM"
Coefficients:
(Intercept) x1b x1c x2B x2C x3 x4 x5
0.4381 0.5968 0.0344 0.2012 0.1789 -0.1320 -0.2155 NA
x1b:x2B x1c:x2B x1b:x2C x1c:x2C
-1.8763 NA -0.8651 NA
> summary(rm1) -> s1
> confint(rm1) -> ci1
> stopifnot(identical(is.na(coef(cm1)), apply(ci1, 1L, anyNA)),
+ identical(sigma(rm1), s1$ sigma),
+ identical(vcov(rm1, complete=FALSE), s1$ cov ),
+ TRUE)
>
> print(s1, showAlgo=FALSE)
Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data)
\--> method = "MM"
Residuals:
Min 1Q Median 3Q Max
-1.4584 -0.3556 0.0246 0.3651 1.0296
Coefficients: (3 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4381 0.5443 0.80 0.44
x1b 0.5968 0.6423 0.93 0.38
x1c 0.0344 0.6880 0.05 0.96
x2B 0.2012 0.7164 0.28 0.79
x2C 0.1789 0.6871 0.26 0.80
x3 -0.1320 0.4155 -0.32 0.76
x4 -0.2155 0.1694 -1.27 0.24
x5 NA NA NA NA
x1b:x2B -1.8763 1.2153 -1.54 0.16
x1c:x2B NA NA NA NA
x1b:x2C -0.8651 0.7466 -1.16 0.28
x1c:x2C NA NA NA NA
Robust residual standard error: 0.927
(3 observations deleted due to missingness)
Multiple R-squared: 0.338, Adjusted R-squared: -0.251
Convergence in 15 IRWLS iterations
Robustness weights:
2 weights are ~= 1. The remaining 16 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.787 0.937 0.985 0.952 0.988 0.994
> ci1
2.5 % 97.5 %
(Intercept) -0.79333 1.66946
x1b -0.85607 2.04973
x1c -1.52188 1.59076
x2B -1.41948 1.82189
x2C -1.37549 1.73320
x3 -1.07182 0.80783
x4 -0.59863 0.16756
x5 NA NA
x1b:x2B -4.62539 0.87283
x1c:x2B NA NA
x1b:x2C -2.55391 0.82381
x1c:x2C NA NA
> ## drop1 should return df = 0
> #drop1(rm1) ## drop.lm does not return valid results (yet)!
>
> ####
> ## methods that should just drop the NA coefficients
> ## m3 is actually the same as rm1, so anova should raise an error
> tools::assertError(anova(rm1, m3, test="Wald"))
> tools::assertError(anova(rm1, m3, test="Deviance"))
> ## but comparing rm1 and rm0 should be ok
> anova(rm1, rm0, test="Wald")
Robust Wald Test Table
Model 1: y ~ x1 * x2 + x3 + x4 + x5
Model 2: y ~ x1 * x2 + x3
Largest model fitted by lmrob(), i.e. SM
pseudoDf Test.Stat Df Pr(>chisq)
1 6
2 10 1.62 1 0.2
> anova(rm1, rm0, test="Deviance")
Robust Deviance Table
Model 1: y ~ x1 * x2 + x3 + x4 + x5
Model 2: y ~ x1 * x2 + x3
Largest model fitted by lmrob(), i.e. SM
pseudoDf Test.Stat Df Pr(>chisq)
1 6
2 10 1.4 1 0.24
> ## commands with single #:
> ## they do (or might) not return sensible results for robust fits
> ## and need to be checked again
> ## IGNORE_RDIFF_BEGIN
> cooks.distance(rm1)
4 5 7 8 10 11 12 13 14 16
0.39378416 0.14209881 0.09280251 0.06910325 0.11270605 0.04346678 0.04339510 0.11965512 0.35871198 0.00060833
17 19 20 21 22 23 25 26
0.02028992 0.11270605 0.04346678 0.04339510 0.00218446 0.02376421 0.01556261 0.01127912
> deviance(rm1)
[1] 7.507
> dfbeta(rm1)
(Intercept) x1b x1c x2B x2C x3 x4 x1b:x2B x1b:x2C
4 0.5298595 -0.613726 -0.679075 0.17409248 -0.637481 0.442340 -0.1267794 -0.1777503 0.702152
5 -0.1947849 0.220978 0.307212 0.19557269 0.102180 -0.038671 0.1222857 0.4270748 -0.182749
7 -0.2193956 0.246729 0.372951 0.20853644 -0.284341 0.014404 0.1731280 -0.1327858 0.167032
8 -0.1566493 0.178390 0.238672 0.16094447 0.101442 -0.049168 0.0873118 -0.1291981 0.219182
10 -0.4229180 0.395149 0.354080 0.40584962 0.424543 0.117095 -0.0663878 -0.4208686 -0.385308
11 -0.0440382 0.438017 0.039727 0.05718669 0.091351 -0.072743 -0.0114323 -0.4560648 -0.465901
12 0.0397872 -0.047201 -0.428171 -0.05112319 -0.079674 0.063041 0.0086918 0.0620558 0.070863
13 0.1412832 -0.161587 -0.206625 -0.69496160 -0.111315 0.062934 -0.0673962 0.6733146 0.153768
14 -0.2188246 0.261364 0.182329 0.29073496 0.488521 -0.393908 -0.0766202 -1.1690942 -0.418435
16 0.0123913 -0.013827 -0.022407 -0.01119215 -0.031096 -0.003704 -0.0115431 0.0058277 0.039046
17 -0.0072068 0.006731 0.029304 -0.00058975 -0.037398 0.037184 0.0281032 0.0168727 -0.197699
19 0.7910973 -0.818866 -0.859935 -0.80816564 -0.789472 0.117095 -0.0663878 0.7931467 0.828707
20 -0.0440382 -0.333308 0.039727 0.05718669 0.091351 -0.072743 -0.0114323 0.3152602 0.305424
21 0.0397872 -0.047201 0.353897 -0.05112319 -0.079674 0.063041 0.0086918 0.0620558 0.070863
22 -0.0257286 0.030017 0.030295 0.09637704 0.037104 -0.027246 0.0026346 -0.0983665 -0.037602
23 0.0180167 -0.017759 -0.061691 -0.00357228 0.066938 -0.068056 -0.0550510 0.1788593 -0.026660
25 -0.0472868 0.053847 0.072079 0.04856922 0.220617 -0.014772 0.0263989 -0.0389600 -0.237709
26 -0.0564761 0.065372 0.072918 0.06375284 0.066714 -0.045991 0.0142193 -0.0629283 -0.229144
> dfbetas(rm1)
(Intercept) x1b x1c x2B x2C x3 x4 x1b:x2B x1b:x2C
4 0.800578 -0.700136 -0.720945 0.21351154 -0.789273 1.4246062 -0.5842907 -0.1668632 0.650705
5 -0.277301 0.237525 0.307309 0.22599707 0.119201 -0.1173486 0.5310176 0.3777527 -0.159574
7 -0.279327 0.237176 0.333641 0.21550902 -0.296648 0.0390900 0.6723410 -0.1050375 0.130435
8 -0.206807 0.177817 0.221401 0.17246923 0.109742 -0.1383605 0.3515988 -0.1059744 0.177481
10 -0.550218 0.388155 0.323685 0.42859027 0.452603 0.3247238 -0.2634535 -0.3401988 -0.307466
11 -0.055794 0.419000 0.035366 0.05880998 0.094839 -0.1964475 -0.0441802 -0.3589978 -0.362044
12 0.050473 -0.045209 -0.381658 -0.05264183 -0.082822 0.1704634 0.0336327 0.0489107 0.055137
13 0.204112 -0.176259 -0.209751 -0.81496353 -0.131780 0.1938028 -0.2969969 0.6043721 0.136256
14 -0.373966 0.337246 0.218944 0.40330287 0.684126 -1.4349146 -0.3994073 -1.2413453 -0.438605
16 0.015390 -0.012966 -0.019555 -0.01128348 -0.031649 -0.0098062 -0.0437311 0.0044971 0.029746
17 -0.009121 0.006432 0.026060 -0.00060585 -0.038785 0.1003120 0.1084907 0.0132676 -0.153466
19 1.029221 -0.804371 -0.786115 -0.85344894 -0.841651 0.3247238 -0.2634535 0.6411206 0.661286
20 -0.055794 -0.318837 0.035366 0.05880998 0.094839 -0.1964475 -0.0441802 0.2481615 0.237340
21 0.050473 -0.045209 0.315452 -0.05264183 -0.082822 0.1704634 0.0336327 0.0489107 0.055137
22 -0.031987 0.028176 0.026465 0.09725830 0.037800 -0.0722026 0.0099908 -0.0759818 -0.028673
23 0.022591 -0.016813 -0.054352 -0.00363580 0.068777 -0.1818947 -0.2105514 0.1393397 -0.020503
25 -0.059657 0.051292 0.063896 0.04973754 0.228076 -0.0397252 0.1015891 -0.0305387 -0.183941
26 -0.070783 0.061862 0.064216 0.06485799 0.068517 -0.1228668 0.0543602 -0.0490026 -0.176150
> if(FALSE)
+ effects(rm1) ## fails
> ## IGNORE_RDIFF_END
> extractAIC(rm1)
[1] 9.0000 2.2583
> infl.1 <- influence(rm1)
> ## checking robustbase:::.lmrob.hat() which uses qr(.)
> hatvals_lm <- stats:::hatvalues.lm # just to check that it's the same computations
> stopifnot(identical(infl.1, lm.influence(rm1)),
+ setequal(names(infl.1), c("hat", "coefficients", "sigma", "wt.res")),
+ all.equal(hv1 <- hatvalues(rm1), .lmrob.hat(wqr=rm1$qr), tol=1e-15),
+ all.equal(hv1, hatvals_lm(rm1, infl.1), tol=1e-15),
+ all.equal(hat(cm1$qr), unname(hatvalues(cm1)), tol=1e-15),
+ all.equal(unname(hv1), hat(rm1$qr), tol=1e-15),
+ ## ditto :
+ all.equal(hv1c <- hatvalues(rm1c), hatvals_lm(rm1c), tol=1e-15))
>
> ## kappa() & labels() :
> stopifnot(is.infinite(kr1 <- kappa(rm1)), kr1 == kappa(cm1), # = +Inf both
+ identical(labels(rm1), labels(cm1)))
> logLik(rm1)# well, and what does it mean?
'log Lik.' -17.67 (df=10)
> ## plot(rm1, which=1) ## plot.lmrob() fails "singular covariance" .. FIXME!
> par(mfrow=c(2,2))
> plot(rm1, which=2:4)
> stopifnot(all.equal(predict(rm1), predict(rm1c), tol=1e-15),
+ all.equal(predict(rm1, se.fit=TRUE, interval="confidence"),
+ predict(rm1c, se.fit=TRUE, interval="confidence"), tol=4e-15)) # seen 1.3e-15 (ATLAS)
> predict(rm1, type="terms", se.fit=TRUE, interval="confidence")
$fit
x1 x2 x3 x4 x5 x1:x2
4 -0.26908 0.074520 -0.166290 0.17233795 0 0.45689
5 0.32774 0.074520 0.026620 -0.03309916 0 -1.41939
7 -0.26908 0.052168 -0.038119 0.28384254 0 0.45689
8 0.32774 0.052168 0.020155 -0.26844357 0 -0.40816
10 -0.26908 -0.126688 0.194821 -0.38642275 0 0.45689
11 0.32774 -0.126688 0.067831 0.11957373 0 0.45689
12 -0.23465 -0.126688 0.065098 0.26547275 0 0.45689
13 -0.26908 0.074520 0.020882 -0.08237063 0 0.45689
14 0.32774 0.074520 -0.132148 0.06953345 0 -1.41939
16 -0.26908 0.052168 -0.087685 -0.47721028 0 0.45689
17 0.32774 0.052168 0.034769 0.04888197 0 -0.40816
19 -0.26908 -0.126688 0.046496 -0.10823918 0 0.45689
20 0.32774 -0.126688 -0.078945 0.03438888 0 0.45689
21 -0.23465 -0.126688 -0.060426 0.20062634 0 0.45689
22 -0.26908 0.074520 0.103967 -0.00026715 0 0.45689
23 0.32774 0.074520 0.106440 0.42945756 0 -1.41939
25 -0.26908 0.052168 -0.035072 -0.27545520 0 0.45689
26 0.32774 0.052168 -0.088392 0.00739277 0 -0.40816
attr(,"constant")
[1] 0.32347
$se.fit
x1 x2 x3 x4 x5 x1:x2
4 0.35192 0.42010 0.523390 0.13540939 0 0.29013
5 0.29582 0.42010 0.083786 0.02600668 0 0.95012
7 0.35192 0.40345 0.119979 0.22302078 0 0.29013
8 0.29582 0.40345 0.063436 0.21092151 0 0.53827
10 0.35192 0.40191 0.613190 0.30362011 0 0.29013
11 0.29582 0.40191 0.213494 0.09395148 0 0.29013
12 0.40411 0.40191 0.204892 0.20858727 0 0.29013
13 0.35192 0.42010 0.065724 0.06472026 0 0.29013
14 0.29582 0.42010 0.415930 0.05463383 0 0.95012
16 0.35192 0.40345 0.275984 0.37495370 0 0.29013
17 0.29582 0.40345 0.109434 0.03840755 0 0.53827
19 0.35192 0.40191 0.146343 0.08504570 0 0.29013
20 0.29582 0.40191 0.248476 0.02702003 0 0.29013
21 0.40411 0.40191 0.190187 0.15763614 0 0.29013
22 0.35192 0.42010 0.327230 0.00020991 0 0.29013
23 0.29582 0.42010 0.335015 0.33743343 0 0.95012
25 0.35192 0.40345 0.110386 0.21643068 0 0.29013
26 0.29582 0.40345 0.278210 0.00580864 0 0.53827
$lwr
x1 x2 x3 x4 x5 x1:x2
4 -1.06517 -0.87582 -1.35028 -0.1339794 0 -0.19943
5 -0.34144 -0.87582 -0.16292 -0.0919303 0 -3.56872
7 -1.06517 -0.86049 -0.30953 -0.2206655 0 -0.19943
8 -0.34144 -0.86049 -0.12335 -0.7455812 0 -1.62581
10 -1.06517 -1.03588 -1.19231 -1.0732591 0 -0.19943
11 -0.34144 -1.03588 -0.41513 -0.0929593 0 -0.19943
12 -1.14880 -1.03588 -0.39840 -0.2063844 0 -0.19943
13 -1.06517 -0.87582 -0.12780 -0.2287780 0 -0.19943
14 -0.34144 -0.87582 -1.07305 -0.0540569 0 -3.56872
16 -1.06517 -0.86049 -0.71200 -1.3254145 0 -0.19943
17 -0.34144 -0.86049 -0.21279 -0.0380019 0 -1.62581
19 -1.06517 -1.03588 -0.28455 -0.3006259 0 -0.19943
20 -0.34144 -1.03588 -0.64104 -0.0267347 0 -0.19943
21 -1.14880 -1.03588 -0.49066 -0.1559714 0 -0.19943
22 -1.06517 -0.87582 -0.63628 -0.0007420 0 -0.19943
23 -0.34144 -0.87582 -0.65142 -0.3338699 0 -3.56872
25 -1.06517 -0.86049 -0.28478 -0.7650554 0 -0.19943
26 -0.34144 -0.86049 -0.71775 -0.0057473 0 -1.62581
attr(,"constant")
[1] 0.32347
$upr
x1 x2 x3 x4 x5 x1:x2
4 0.52701 1.02486 1.01770 0.47865527 0 1.11321
5 0.99693 1.02486 0.21616 0.02573203 0 0.72993
7 0.52701 0.96483 0.23329 0.78835059 0 1.11321
8 0.99693 0.96483 0.16366 0.20869402 0 0.80949
10 0.52701 0.78250 1.58195 0.30041366 0 1.11321
11 0.99693 0.78250 0.55079 0.33210673 0 1.11321
12 0.67950 0.78250 0.52860 0.73732993 0 1.11321
13 0.52701 1.02486 0.16956 0.06403677 0 1.11321
14 0.99693 1.02486 0.80875 0.19312376 0 0.72993
16 0.52701 0.96483 0.53663 0.37099391 0 1.11321
17 0.99693 0.96483 0.28233 0.13576588 0 0.80949
19 0.52701 0.78250 0.37755 0.08414755 0 1.11321
20 0.99693 0.78250 0.48315 0.09551244 0 1.11321
21 0.67950 0.78250 0.36981 0.55722407 0 1.11321
22 0.52701 1.02486 0.84421 0.00020769 0 1.11321
23 0.99693 1.02486 0.86430 1.19278501 0 0.72993
25 0.52701 0.96483 0.21464 0.21414502 0 1.11321
26 0.99693 0.96483 0.54096 0.02053283 0 0.80949
attr(,"constant")
[1] 0.32347
$df
[1] 9
$residual.scale
[1] 0.92726
> ## proj(rm1) ## --> effects() [see FIXME above]: fails
> residuals(rm1)
4 5 7 8 10 11 12 13 14 16 17
1.003436 1.029645 -0.321738 0.691394 -0.498376 0.342960 -0.359752 -1.145548 -1.458427 -0.043483 -0.395061
19 20 21 22 23 25 26
0.498376 -0.342960 0.359752 0.092640 0.232325 0.366908 -0.270349
> ## the next two work via lm.influence() == infl.1
> ## IGNORE_RDIFF_BEGIN
> stopifnot(identical(print(rstandard(rm1)), rstandard(rm1, infl.1)),
+ identical(print(rstudent (rm1)), rstudent (rm1, infl.1)))
4 5 7 8 10 11 12 13 14 16 17
1.603821 1.376940 -0.622406 0.956366 -0.840478 0.559647 -0.576799 -1.438932 -1.934023 -0.069523 -0.545099
19 20 21 22 23 25 26
0.840478 -0.559647 0.576799 0.142326 0.392184 0.498702 -0.383397
4 5 7 8 10 11 12 13 14 16 17
1.961168 1.586452 -0.641320 1.021827 -0.884954 0.573836 -0.592182 -1.682423 -2.674934 -0.069884 -0.558330
19 20 21 22 23 25 26
0.884954 -0.573836 0.592182 0.143204 0.397981 0.509193 -0.388893
> ## IGNORE_RDIFF_END
> simulate(rm1)
sim_1
4 0.494959
5 -0.604246
7 2.651906
8 0.793257
10 0.976483
11 1.884235
12 2.556984
13 1.100485
14 -0.029195
16 0.802217
17 1.192705
19 -0.250731
20 0.178226
21 1.827698
22 1.309304
23 1.482225
25 -0.996092
26 2.375265
> V1 <- vcov(rm1, complete=FALSE)
> ## but don't show the "eigen" part {vectors may flip sign}:
> attributes(V1) <- attributes(V1)[c("dim","dimnames", "weights")]; V1
(Intercept) x1b x1c x2B x2C x3 x4 x1b:x2B x1b:x2C
(Intercept) 0.296312 -0.321429 -0.338842 -0.238010 -0.3289125 0.1357438 -0.0352579 0.305162 0.321423
x1b -0.321429 0.412501 0.369763 0.253038 0.3616767 -0.1594475 0.0395980 -0.399754 -0.414159
x1c -0.338842 0.369763 0.473317 0.274811 0.3497592 -0.1464335 0.0605871 -0.273260 -0.349092
x2B -0.238010 0.253038 0.274811 0.513277 0.2342086 -0.0640599 0.0353593 -0.557840 -0.233097
x2C -0.328913 0.361677 0.349759 0.234209 0.4721185 -0.2294044 0.0024864 -0.521954 -0.439408
x3 0.135744 -0.159448 -0.146434 -0.064060 -0.2294044 0.1726038 0.0037187 0.308735 0.203925
x4 -0.035258 0.039598 0.060587 0.035359 0.0024864 0.0037187 0.0286797 0.063743 -0.012435
x1b:x2B 0.305162 -0.399754 -0.273260 -0.557840 -0.5219539 0.3087350 0.0637434 1.476860 0.498060
x1b:x2C 0.321423 -0.414159 -0.349092 -0.233097 -0.4394078 0.2039253 -0.0124347 0.498060 0.557368
attr(,"weights")
4 5 7 8 10 11 12 13 14 16 17 19 20
0.89614 0.89081 0.98906 0.94998 0.97385 0.98757 0.98633 0.86577 0.78729 0.99980 0.98353 0.97385 0.98757
21 22 23 25 26
0.98633 0.99909 0.99429 0.98578 0.99227
> set.seed(12); sc <- simulate(cm1, 64)
> set.seed(12); rc <- simulate(rm1, 64)
>
> stopifnot(all.equal(sqrt(diag(V1)), coef(summary(rm1))[,"Std. Error"], tol=1e-15),
+ all.equal(sc, rc, tolerance = 0.08),# dimension *and* approx. values (no NA)
+ identical(variable.names(rm1), variable.names(cm1)),
+ all.equal(residuals(rm1), residuals(cm1), tolerance = 0.05),# incl. names
+ all.equal(rstudent (rm1), rstudent (cm1), tolerance = 0.06),
+ identical(dimnames(rm1), dimnames(cm1)),
+ all.equal(dummy.coef(rm1), dummy.coef(cm1), tolerance= .5)) ## check mostly structure
>
> ## other helper functions
> stopifnot(identical(case.names(rm1), case.names(cm1)),
+ all.equal(family(rm1), family(cm1)),# identical() upto environment
+ identical(formula(rm1), formula(cm1)),
+ nobs(rm1) == nobs(cm1))
> #add1(rm0, ~ . + x3 + x4 + x5) ## does not return valid results (yet)!
>
>
> ## test other initial estimators
> lmrob(y ~ x1*x2 + x3 + x4 + x5, data, init="M-S")
Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data, init = "M-S")
\--> method = "M-SM"
Coefficients:
(Intercept) x1b x1c x2B x2C x3 x4 x5
0.4358 0.5996 0.0346 0.2005 0.1877 -0.1395 -0.2185 NA
x1b:x2B x1c:x2B x1b:x2C x1c:x2C
-1.8957 NA -0.8698 NA
Warning message:
In lmrob.M.S(x, y, control, mf = mf) :
Skipping design matrix equilibration (DGEEQU): row 12 is exactly zero.
> lmrob(y ~ x1*x2 + x3 + x4 + x5, data, init=lmrob.lar)
Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data, init = lmrob.lar)
\--> method = "lM"
Coefficients:
(Intercept) x1b x1c x2B x2C x3 x4 x5
0.561131 0.444339 0.000184 0.530303 -0.251794 0.236541 -0.082680 NA
x1b:x2B x1c:x2B x1b:x2C x1c:x2C
-1.298418 NA -0.597602 NA
>
> ## test all zero design matrix
> data <- data.frame(y=1:10,x1=0,x2=0,os=2,w=c(0.5, 1))
> (m5 <- lmrob(y ~ 1+x1+x2+offset(os), data, weights=w))
Call:
lmrob(formula = y ~ 1 + x1 + x2 + offset(os), data = data, weights = w)
\--> method = "MM"
Coefficients:
(Intercept) x1 x2
3.64 NA NA
> (sm5 <- summary(m5))
Call:
lmrob(formula = y ~ 1 + x1 + x2 + offset(os), data = data, weights = w)
\--> method = "MM"
Weighted Residuals:
Min 1Q Median 3Q Max
-3.6412 -1.8110 -0.0473 2.0093 4.3588
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.64 1.03 3.53 0.0064 **
x1 NA NA NA NA
x2 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Robust residual standard error: 3.24
Convergence in 8 IRWLS iterations
Robustness weights:
1 2 3 4 5 6 7 8 9 10
0.909 0.889 0.970 0.977 0.998 0.999 0.992 0.952 0.952 0.842
Algorithmic parameters:
tuning.chi bb tuning.psi refine.tol rel.tol scale.tol
1.55e+00 5.00e-01 4.69e+00 1.00e-07 1.00e-07 1.00e-10
solve.tol zero.tol eps.outlier eps.x warn.limit.reject warn.limit.meanrw
1.00e-07 1.00e-10 1.00e-02 1.82e-12 5.00e-01 5.00e-01
nResample max.it best.r.s k.fast.s k.max maxit.scale trace.lev
500 50 2 1 200 200 0
mts compute.rd fast.s.large.n
1000 0 2000
psi subsampling cov compute.outlier.stats
"bisquare" "nonsingular" ".vcov.avar1" "SM"
seed : int(0)
> (m6 <- lmrob(y ~ 0+x1+x2+offset(os), data, weights=w))
Call:
lmrob(formula = y ~ 0 + x1 + x2 + offset(os), data = data, weights = w)
\--> method = "MM"
Coefficients:
x1 x2
NA NA
> (sm6 <- summary(m6))
Call:
lmrob(formula = y ~ 0 + x1 + x2 + offset(os), data = data, weights = w)
\--> method = "MM"
Weighted Residuals:
Min 1Q Median 3Q Max
-3.83 -1.37 1.09 3.54 6.00
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
x1 NA NA NA NA
x2 NA NA NA NA
Robust residual standard error: NA
Convergence in 0 IRWLS iterations
Robustness weights:
[1] NA NA NA NA NA NA NA NA NA NA
Algorithmic parameters:
tuning.psi rel.tol scale.tol solve.tol zero.tol eps.outlier
4.69e+00 1.00e-07 1.00e-10 1.00e-07 1.00e-10 1.00e-02
warn.limit.reject warn.limit.meanrw
5.00e-01 5.00e-01
max.it maxit.scale trace.lev compute.rd fast.s.large.n eps.x
50 200 0 0 2000 0
psi cov compute.outlier.stats
"bisquare" ".vcov.avar1" "SM"
seed : int(0)
>
> sc5 <- summary(cm5 <- lm(y ~ 1+x1+x2+offset(os), data, weights=w))
> sc6 <- summary(cm6 <- lm(y ~ 0+x1+x2+offset(os), data, weights=w))
>
> if(getRversion() <= "3.5.1" && as.numeric(R.version$`svn rev`) < 74993)
+ ## in the past, lm() returned logical empty matrix
+ storage.mode(sc6$coefficients) <- "double"
>
> stopifnot(all.equal(coef(m5), coef(cm5), tolerance = 0.01),
+ all.equal(coef(m6), coef(cm6), tolerance = 1e-14),
+ all.equal(coef(sm5), coef(sc5), tolerance = 0.05),
+ all.equal(coef(sm6), coef(sc6), tolerance = 1e-14),
+ identical(sm5$df, sc5$df),
+ identical(sm6$df, sc6$df))
>
> options(op) # revert
>
> proc.time()
user system elapsed
0.437 0.311 1.059
|
