R Under development (unstable) (2026-02-03 r89375 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. 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.31 0.06 0.37