Running ‘LTS-specials.R’
  Running ‘MCD-specials.R’
  Comparing ‘MCD-specials.Rout’ to ‘MCD-specials.Rout.save’ ... OK
  Running ‘MT-tst.R’ [13s/14s]
  Running ‘NAcoef.R’
  Running ‘OGK-ex.R’
  Comparing ‘OGK-ex.Rout’ to ‘OGK-ex.Rout.save’ ... OK
  Running ‘Rsquared.R’
  Running ‘binom-ni-small.R’
  Comparing ‘binom-ni-small.Rout’ to ‘binom-ni-small.Rout.save’ ... OK
  Running ‘binom-no-x.R’
  Running ‘comedian-tst.R’
  Running ‘exact-fit-categorical.R’
  Running ‘glmrob-1.R’ [15s/15s]
  Running ‘glmrob-specials.R’
  Running ‘huber-etc.R’
  Comparing ‘huber-etc.Rout’ to ‘huber-etc.Rout.save’ ... OK
  Running ‘lmrob-data.R’ [24s/24s]
  Running ‘lmrob-ex12.R’
  Running ‘lmrob-methods.R’
  Comparing ‘lmrob-methods.Rout’ to ‘lmrob-methods.Rout.save’ ... OK
  Running ‘lmrob-outlierStats.R’
  Running ‘lmrob-psifns.R’
  Comparing ‘lmrob-psifns.Rout’ to ‘lmrob-psifns.Rout.save’ ... OK
  Running ‘m-s-estimator.R’
  Running ‘mc-etc.R’
  Running ‘mc-strict.R’ [14s/14s]
  Running ‘nlregrob-tst.R’ [43s/44s]
  Running ‘nlrob-tst.R’
  Running ‘plot-ex.R’
  Running ‘poisson-ex.R’
  Running ‘psi-rho-etc.R’
  Comparing ‘psi-rho-etc.Rout’ to ‘psi-rho-etc.Rout.save’ ... OK
  Running ‘small-sample.R’
  Comparing ‘small-sample.Rout’ to ‘small-sample.Rout.save’ ... OK
  Running ‘subsample.R’
  Running ‘tlts.R’
  Comparing ‘tlts.Rout’ to ‘tlts.Rout.save’ ... OK
  Running ‘tmcd.R’ [12s/12s]
  Running ‘weights.R’
  Comparing ‘weights.Rout’ to ‘weights.Rout.save’ ... OK
  Running ‘wgt-himed-xtra.R’
  Running ‘wgt-himed.R’
  Comparing ‘wgt-himed.Rout’ to ‘wgt-himed.Rout.save’ ... OK
Running the tests in ‘tests/NAcoef.R’ failed.
Complete output:
  > ## test handing 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
  > source(system.file("test-tools-1.R", package="Matrix", mustWork=TRUE))
  Loading required package: tools
  > ##-> assertError(), etc
  > options(digits = 5)# -> 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  
       0.4381       0.5968       0.0344       0.2012       0.1789      -0.1320  
           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
      -0.2155           NA      -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 
  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
  > assertError(anova(rm1, m3, test="Wald"))
  > 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
  > #cooks.distance(rm1)
  > #deviance(rm1)
  > #dfbeta(rm1)
  > #dfbetas(rm1)
  > #effects(rm1) ## fails
  > #extractAIC(rm1)
  > #influence(rm1)
  > stopifnot(all.equal(hv1 <- hatvalues(rm1), .lmrob.hat(wqr=rm1$qr), tol=1e-15),
  +           all.equal(hv1, stats:::hatvalues.lm(rm1), 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), stats:::hatvalues.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=1e-15))
  > 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) ## fails "FIXME"
  > residuals(rm1)
          4         5         7         8        10        11        12        13 
   1.003436  1.029645 -0.321738  0.691394 -0.498376  0.342960 -0.359752 -1.145548 
         14        16        17        19        20        21        22        23 
  -1.458427 -0.043483 -0.395061  0.498376 -0.342960  0.359752  0.092640  0.232325 
         25        26 
   0.366908 -0.270349 
  > #rstandard(rm1)
  > #rstudent(rm1)
  > #simulate(rm1) ## just $weights needs to be changed to prior weights
  > 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
  (Intercept)    0.296312 -0.321429 -0.338842 -0.238010 -0.3289125  0.1357438
  x1b           -0.321429  0.412501  0.369763  0.253038  0.3616767 -0.1594475
  x1c           -0.338842  0.369763  0.473317  0.274811  0.3497592 -0.1464335
  x2B           -0.238010  0.253038  0.274811  0.513277  0.2342086 -0.0640599
  x2C           -0.328913  0.361677  0.349759  0.234209  0.4721185 -0.2294044
  x3             0.135744 -0.159448 -0.146434 -0.064060 -0.2294044  0.1726038
  x4            -0.035258  0.039598  0.060587  0.035359  0.0024864  0.0037187
  x1b:x2B        0.305162 -0.399754 -0.273260 -0.557840 -0.5219539  0.3087350
  x1b:x2C        0.321423 -0.414159 -0.349092 -0.233097 -0.4394078  0.2039253
                      x4   x1b:x2B   x1b:x2C
  (Intercept) -0.0352579  0.305162  0.321423
  x1b          0.0395980 -0.399754 -0.414159
  x1c          0.0605871 -0.273260 -0.349092
  x2B          0.0353593 -0.557840 -0.233097
  x2C          0.0024864 -0.521954 -0.439408
  x3           0.0037187  0.308735  0.203925
  x4           0.0286797  0.063743 -0.012435
  x1b:x2B      0.0637434  1.476860  0.498060
  x1b:x2C     -0.0124347  0.498060  0.557368
  attr(,"weights")
        4       5       7       8      10      11      12      13      14      16 
  0.89614 0.89081 0.98906 0.94998 0.97385 0.98757 0.98633 0.86577 0.78729 0.99980 
       17      19      20      21      22      23      25      26 
  0.98353 0.97385 0.98757 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  
       0.4358       0.5996       0.0346       0.2005       0.1877      -0.1395  
           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
      -0.2185           NA      -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  
     0.561131     0.444339     0.000184     0.530303    -0.251794     0.236541  
           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
    -0.082680           NA    -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"
  Residuals:
     Min     1Q Median     3Q    Max 
  -4.641 -2.391 -0.141  2.109  4.359 
  
  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 
           1.55e+00          5.00e-01          4.69e+00          1.00e-07 
            rel.tol         scale.tol         solve.tol       eps.outlier 
           1.00e-07          1.00e-10          1.00e-07          1.00e-02 
              eps.x warn.limit.reject warn.limit.meanrw 
           1.82e-12          5.00e-01          5.00e-01 
       nResample         max.it       best.r.s       k.fast.s          k.max 
             500             50              2              1            200 
     maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
             200              0           1000              0           2000 
                    psi           subsampling                   cov 
             "bisquare"         "nonsingular"         ".vcov.avar1" 
  compute.outlier.stats 
                   "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 
           4.69e+00          1.00e-07          1.00e-10          1.00e-07 
        eps.outlier warn.limit.reject warn.limit.meanrw 
           1.00e-02          5.00e-01          5.00e-01 
          max.it    maxit.scale      trace.lev     compute.rd fast.s.large.n 
              50            200              0              0           2000 
           eps.x 
               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))
  > 
  > stopifnot(all.equal(coef(m5), coef(cm5), tolerance = 0.01),
  +           identical(coef(m6), coef(cm6)),
  +           all.equal(coef(sm5), coef(sc5), tolerance = 0.05),
  +           identical(coef(sm6), coef(sc6)),
  +           identical(sm5$df, sc5$df),
  +           identical(sm6$df, sc6$df))
  Error: identical(coef(sm6), coef(sc6)) is not TRUE
  Execution halted
Running the tests in ‘tests/Rsquared.R’ failed.
Complete output:
  > require(robustbase)
  Loading required package: robustbase
  > 
  > set.seed(17)# reproducibility!
  > ## to check:
  > ## - for the empty model
  > summary(lmrob(Y ~ 0, coleman))
  
  Call:
  lmrob(formula = Y ~ 0, data = coleman)
   \--> method = ""
  Residuals:
     Min     1Q Median     3Q    Max 
   22.70  32.78  35.86  39.95  43.10 
  
  No Coefficients
  > ## - with and without an intercept in the  model
  > summary(lmrob(Y ~ 1, coleman))
  
  Call:
  lmrob(formula = Y ~ 1, data = coleman)
   \--> method = "MM"
  Residuals:
       Min       1Q   Median       3Q      Max 
  -12.8605  -2.7855   0.2945   4.3895   7.5395 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   35.560      1.342    26.5   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Robust residual standard error: 5.48 
  Convergence in 9 IRWLS iterations
  
  Robustness weights: 
   one weight is ~= 1. The remaining 19 ones are summarized as
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.5611  0.8986  0.9553  0.9044  0.9918  0.9987 
  Algorithmic parameters: 
         tuning.chi                bb        tuning.psi        refine.tol 
          1.548e+00         5.000e-01         4.685e+00         1.000e-07 
            rel.tol         scale.tol         solve.tol       eps.outlier 
          1.000e-07         1.000e-10         1.000e-07         5.000e-03 
              eps.x warn.limit.reject warn.limit.meanrw 
          1.819e-12         5.000e-01         5.000e-01 
       nResample         max.it       best.r.s       k.fast.s          k.max 
             500             50              2              1            200 
     maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
             200              0           1000              0           2000 
                    psi           subsampling                   cov 
             "bisquare"         "nonsingular"         ".vcov.avar1" 
  compute.outlier.stats 
                   "SM" 
  seed : int(0) 
  > writeLines(sfm <- capture.output(
  +                        summary(lmrob(Y ~ ., coleman)))) # and this must be "identical":
  
  Call:
  lmrob(formula = Y ~ ., data = coleman)
   \--> method = "MM"
  Residuals:
       Min       1Q   Median       3Q      Max 
  -4.16181 -0.39226  0.01611  0.55619  7.22766 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 30.50232    6.71260   4.544 0.000459 ***
  salaryP     -1.66615    0.43129  -3.863 0.001722 ** 
  fatherWc     0.08425    0.01467   5.741 5.10e-05 ***
  sstatus      0.66774    0.03385  19.726 1.30e-11 ***
  teacherSc    1.16778    0.10983  10.632 4.35e-08 ***
  motherLev   -4.13657    0.92084  -4.492 0.000507 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Robust residual standard error: 1.134 
  Multiple R-squared:  0.9814,	Adjusted R-squared:  0.9747 
  Convergence in 11 IRWLS iterations
  
  Robustness weights: 
   observation 18 is an outlier with |weight| = 0 ( < 0.005); 
   The remaining 19 ones are summarized as
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.1491  0.9412  0.9847  0.9279  0.9947  0.9982 
  Algorithmic parameters: 
         tuning.chi                bb        tuning.psi        refine.tol 
          1.548e+00         5.000e-01         4.685e+00         1.000e-07 
            rel.tol         scale.tol         solve.tol       eps.outlier 
          1.000e-07         1.000e-10         1.000e-07         5.000e-03 
              eps.x warn.limit.reject warn.limit.meanrw 
          1.569e-10         5.000e-01         5.000e-01 
       nResample         max.it       best.r.s       k.fast.s          k.max 
             500             50              2              1            200 
     maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
             200              0           1000              0           2000 
                    psi           subsampling                   cov 
             "bisquare"         "nonsingular"         ".vcov.avar1" 
  compute.outlier.stats 
                   "SM" 
  seed : int(0) 
  > 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))
  
  Call:
  lmrob(formula = Y ~ . - 1, data = coleman, model = FALSE, x = FALSE, y = FALSE)
   \--> method = "MM"
  Residuals:
       Min       1Q   Median       3Q      Max 
  -4.86146 -0.59195 -0.04679  0.87826  5.40639 
  
  Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
  salaryP   -1.97540    0.45262  -4.364 0.000555 ***
  fatherWc   0.03388    0.02220   1.526 0.147749    
  sstatus    0.55922    0.07590   7.367 2.34e-06 ***
  teacherSc  1.60446    0.19039   8.427 4.51e-07 ***
  motherLev -0.48903    0.90805  -0.539 0.598097    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Robust residual standard error: 1.344 
  Multiple R-squared:  0.9987,	Adjusted R-squared:  0.9983 
  Convergence in 14 IRWLS iterations
  
  Robustness weights: 
   3 weights are ~= 1. The remaining 17 ones are summarized as
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.06899 0.89030 0.95860 0.82750 0.98700 0.99820 
  Algorithmic parameters: 
         tuning.chi                bb        tuning.psi        refine.tol 
          1.548e+00         5.000e-01         4.685e+00         1.000e-07 
            rel.tol         scale.tol         solve.tol       eps.outlier 
          1.000e-07         1.000e-10         1.000e-07         5.000e-03 
              eps.x warn.limit.reject warn.limit.meanrw 
          1.569e-10         5.000e-01         5.000e-01 
       nResample         max.it       best.r.s       k.fast.s          k.max 
             500             50              2              1            200 
     maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
             200              0           1000              0           2000 
                    psi           subsampling                   cov 
             "bisquare"         "nonsingular"         ".vcov.avar1" 
  compute.outlier.stats 
                   "SM" 
  seed : int(0) 
  > 
  > ## - 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))
  
  Call:
  lmrob(formula = Y ~ . - 1, data = coleman, weights = wts, model = FALSE, 
      x = FALSE, y = FALSE)
   \--> method = "MM"
  Residuals:
       Min       1Q   Median       3Q      Max 
  -4.53960 -0.03032  0.10984  1.55271  6.36480 
  
  Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
  salaryP   -2.68644    0.05871 -45.761  < 2e-16 ***
  fatherWc   0.04761    0.00721   6.603 8.39e-06 ***
  sstatus    0.58362    0.00314 185.842  < 2e-16 ***
  teacherSc  1.77115    0.07918  22.369 6.20e-13 ***
  motherLev -1.03171    0.34154  -3.021   0.0086 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Robust residual standard error: 0.423 
  Multiple R-squared:  0.9985,	Adjusted R-squared:  0.998 
  Convergence in 5 IRWLS iterations
  
  Robustness weights: 
   3 observations c(12,16,18) are outliers with |weight| = 0 ( < 0.005); 
   5 weights are ~= 1. The remaining 12 ones are summarized as
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.5441  0.9257  0.9833  0.9303  0.9956  0.9985 
  Algorithmic parameters: 
         tuning.chi                bb        tuning.psi        refine.tol 
          1.548e+00         5.000e-01         4.685e+00         1.000e-07 
            rel.tol         scale.tol         solve.tol       eps.outlier 
          1.000e-07         1.000e-10         1.000e-07         5.000e-03 
              eps.x warn.limit.reject warn.limit.meanrw 
          2.219e-10         5.000e-01         5.000e-01 
       nResample         max.it       best.r.s       k.fast.s          k.max 
             500             50              2              1            200 
     maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
             200              0           1000              0           2000 
                    psi           subsampling                   cov 
             "bisquare"         "nonsingular"         ".vcov.avar1" 
  compute.outlier.stats 
                   "SM" 
  seed : int(0) 
  > ## - 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"
  +     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)
  
  Call:
  lmrob(formula = formula, data = data, control = lmrCtrl)
   \--> method = ""
  Residuals:
     Min     1Q Median     3Q    Max 
   22.70  32.78  35.86  39.95  43.10 
  
  No Coefficients
  $residuals
      1     2     3     4     5     6     7     8     9    10    11    12    13 
  37.01 26.51 36.51 40.70 37.10 33.90 41.80 33.40 41.01 37.20 23.30 35.20 34.90 
     14    15    16    17    18    19    20 
  33.10 22.70 39.70 31.80 31.70 43.10 41.01 
  
  $residuals
      1     2     3     4     5     6     7     8     9    10    11    12    13 
  37.01 26.51 36.51 40.70 37.10 33.90 41.80 33.40 41.01 37.20 23.30 35.20 34.90 
     14    15    16    17    18    19    20 
  33.10 22.70 39.70 31.80 31.70 43.10 41.01 
  
  $df
  [1]  0 20  0
  
  $df
  [1]  0 20  0
  
  $coefficients
       Estimate Std. Error t value Pr(>|t|)
  
  $coefficients
       Estimate Std. Error t value Pr(>|t|)
  
  $r.squared
  [1] 0
  
  $r.squared
  [1] 0
  
  $adj.r.squared
  [1] 0
  
  $adj.r.squared
  [1] 0
  
  [1] "Component \"coefficients\": Modes: numeric, logical"
  Error in test(Y ~ 0, coleman, c("residuals", "df", "coefficients", "r.squared",  : 
    all.equal(sc[items], sr[items], tol.. = 1e-10) are not all TRUE
  Execution halted