R version 4.4.1 Patched (2024-06-15 r86768) -- "Race for Your Life"
Copyright (C) 2024 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.

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 handing of weights and offset argument
> require(robustbase)
Loading required package: robustbase
> 
> ## generate simple example data  (extension of the one in ./NAcoef.R )
> data <- expand.grid(x1=letters[1:3], x2=LETTERS[1:4], rep=1:3,
+                     KEEP.OUT.ATTRS = FALSE)
> ## generate offset column
> data$os <- 1:nrow(data)
> set.seed(1)
> data$y <- data$os + rnorm(nrow(data))
> ## add collinear variables
> data$x3 <- rnorm(nrow(data))
> data$x4 <- rnorm(nrow(data))
> data$x5 <- data$x3 + data$x4 ## lm() will have 'x5'  "aliased" (and give coef = NA)
> ## add some NA terms
> data$y[1] <- NA
> data$x4[2:3] <- NA ## to test anova
> ## generate weights
> ## some obs with weight 0
> data$weights <- as.numeric(with(data, x1 != 'c' | (x2 != 'B' & x2 != 'C')))
> ## some obs with weight 2
> data$weights[data$x1 == 'b'] <- 2
> ## data2 := {data + weights}, encoded in "data2" (-> "ok" for coef(), not for SE)
> data2 <- rbind(subset(data, weights  > 0),
+ 	       subset(data, weights == 2))
> ## using these parameters we're essentially forcing lmrob() to
> ## fit a classic model --> easier to compare to lm()
> ctrl <- lmrob.control(psi="optimal", tuning.chi = 20, bb = 0.0003846154,
+                       tuning.psi=20, method="SM", cov=".vcov.w")
> ## SM = MM == the case where  .vcov.avar1  was also defined for
> 
> ## Classical models start with 'cm', robust just with  'rm' (or just 'm'),
> ## "." - models drop 'x5' (which is aliased / extraneous by construction) :
> (cm0 <- lm   (y ~ x1*x2 + x3 + x4 + x5 + offset(os), data))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
    0.01008     -1.14140      0.48156      0.01357      0.86985      0.15178  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
   -0.01655     -0.02388           NA      1.05416     -0.32889      0.69954  
    x1c:x2C      x1b:x2D      x1c:x2D  
   -0.73949      1.08478     -1.31578  

> (cm0.<- lm   (y ~ x1*x2 + x3 + x4      + offset(os), data))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4 + offset(os), data = data)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
    0.01008     -1.14140      0.48156      0.01357      0.86985      0.15178  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
   -0.01655     -0.02388      1.05416     -0.32889      0.69954     -0.73949  
    x1b:x2D      x1c:x2D  
    1.08478     -1.31578  

> (cm1 <- lm   (y ~ x1*x2 + x3 + x4 + x5 + offset(os), data,  weights=weights))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data, 
    weights = weights)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
  -0.021632    -0.079147           NA     1.040529           NA     0.736944  
    x1c:x2C      x1b:x2D      x1c:x2D  
         NA     1.099090    -1.371953  

> (cm1.<- lm   (y ~ x1*x2 + x3 + x4      + offset(os), data,  weights=weights))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4 + offset(os), data = data, 
    weights = weights)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
  -0.021632    -0.079147     1.040529           NA     0.736944           NA  
    x1b:x2D      x1c:x2D  
   1.099090    -1.371953  

> (cm2 <- lm   (y ~ x1*x2 + x3 + x4 + x5,              data2, offset=os))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data2, offset = os)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
  -0.021632    -0.079147           NA     1.040529           NA     0.736944  
    x1c:x2C      x1b:x2D      x1c:x2D  
         NA     1.099090    -1.371953  

> (cm2.<- lm   (y ~ x1*x2 + x3 + x4,                   data2, offset=os))

Call:
lm(formula = y ~ x1 * x2 + x3 + x4, data = data2, offset = os)

Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
  -0.021632    -0.079147     1.040529           NA     0.736944           NA  
    x1b:x2D      x1c:x2D  
   1.099090    -1.371953  

> (rm0 <- lmrob(y ~ x1*x2 + x3 + x4 + x5 + offset(os), data,                   control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data, control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
    0.01008     -1.14140      0.48156      0.01357      0.86985      0.15178  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
   -0.01655     -0.02388           NA      1.05416     -0.32889      0.69954  
    x1c:x2C      x1b:x2D      x1c:x2D  
   -0.73949      1.08478     -1.31578  

> (rm0.<- lmrob(y ~ x1*x2 + x3 + x4      + offset(os), data,                   control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + offset(os), data = data, control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
    0.01008     -1.14140      0.48156      0.01357      0.86985      0.15178  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
   -0.01655     -0.02388      1.05416     -0.32889      0.69954     -0.73949  
    x1b:x2D      x1c:x2D  
    1.08478     -1.31578  

> set.seed(2)
> (rm1 <- lmrob(y ~ x1*x2 + x3 + x4 + x5 + offset(os), data,  weights=weights, control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data, weights = weights,     control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
  -0.021632    -0.079147           NA     1.040529           NA     0.736944  
    x1c:x2C      x1b:x2D      x1c:x2D  
         NA     1.099090    -1.371953  

> set.seed(2)
> (rm1.<- lmrob(y ~ x1*x2 + x3 + x4      + offset(os), data,  weights=weights, control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + offset(os), data = data, weights = weights,     control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
  -0.021632    -0.079147     1.040529           NA     0.736944           NA  
    x1b:x2D      x1c:x2D  
   1.099090    -1.371953  

> set.seed(2)
> (rm2 <- lmrob(y ~ x1*x2 + x3 + x4 + x5,              data2, offset=os,       control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data2, offset = os, control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4           x5      x1b:x2B      x1c:x2B      x1b:x2C  
  -0.021632    -0.079147           NA     1.040529           NA     0.736944  
    x1c:x2C      x1b:x2D      x1c:x2D  
         NA     1.099090    -1.371953  

> set.seed(2)
> (rm2.<- lmrob(y ~ x1*x2 + x3 + x4,                   data2, offset=os,       control=ctrl))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4, data = data2, offset = os, control = ctrl)
 \--> method = "MM"
Coefficients:
(Intercept)          x1b          x1c          x2B          x2C          x2D  
  -0.002961    -1.132857     0.492904     0.017959     0.858031     0.208510  
         x3           x4      x1b:x2B      x1c:x2B      x1b:x2C      x1c:x2C  
  -0.021632    -0.079147     1.040529           NA     0.736944           NA  
    x1b:x2D      x1c:x2D  
   1.099090    -1.371953  

> 
> sc0 <- summary(cm0)
> sc1 <- summary(cm1)
> sc2 <- summary(cm2)
> (sr0 <- summary(rm0))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data, control = ctrl)
 \--> method = "MM"
Residuals:
     Min       1Q   Median       3Q      Max 
-1.50524 -0.48219  0.01663  0.42714  1.59122 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.01008    0.76421   0.013    0.990
x1b         -1.14140    1.02228  -1.117    0.278
x1c          0.48156    1.01891   0.473    0.642
x2B          0.01357    0.95276   0.014    0.989
x2C          0.86985    0.94762   0.918    0.370
x2D          0.15178    0.99480   0.153    0.880
x3          -0.01655    0.22284  -0.074    0.942
x4          -0.02388    0.25629  -0.093    0.927
x5                NA         NA      NA       NA
x1b:x2B      1.05416    1.30705   0.807    0.430
x1c:x2B     -0.32889    1.30044  -0.253    0.803
x1b:x2C      0.69954    1.37279   0.510    0.616
x1c:x2C     -0.73949    1.30141  -0.568    0.577
x1b:x2D      1.08478    1.32102   0.821    0.422
x1c:x2D     -1.31578    1.33335  -0.987    0.336

Robust residual standard error: 1.007 
  (3 observations deleted due to missingness)
Multiple R-squared:  0.9933,	Adjusted R-squared:  0.9887 
Convergence in 1 IRWLS iterations

Robustness weights: 
 All 33 weights are ~= 1.
Algorithmic parameters: 
               bb        refine.tol           rel.tol         scale.tol 
        3.846e-04         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         3.030e-03         4.369e-12 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample     tuning.chi     tuning.psi         max.it       best.r.s 
           500             20             20             50              2 
      k.fast.s          k.max    maxit.scale      trace.lev            mts 
             1            200            200              0           1000 
    compute.rd fast.s.large.n 
             0           2000 
                  psi           subsampling                   cov 
            "optimal"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
> (sr1 <- summary(rm1))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5 + offset(os), data = data, weights = weights, 
    control = ctrl)
 \--> method = "MM"
Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.0956 -0.5369  0.0000  0.3925  2.0381 

Coefficients: (3 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002961   0.977109  -0.003    0.998
x1b         -1.132857   1.133342  -1.000    0.333
x1c          0.492904   1.297399   0.380    0.709
x2B          0.017959   1.213927   0.015    0.988
x2C          0.858031   1.204169   0.713    0.487
x2D          0.208510   1.275792   0.163    0.872
x3          -0.021632   0.284226  -0.076    0.940
x4          -0.079147   0.324629  -0.244    0.811
x5                 NA         NA      NA       NA
x1b:x2B      1.040529   1.443384   0.721    0.482
x1c:x2B            NA         NA      NA       NA
x1b:x2C      0.736944   1.530596   0.481    0.637
x1c:x2C            NA         NA      NA       NA
x1b:x2D      1.099090   1.461384   0.752    0.464
x1c:x2D     -1.371953   1.698858  -0.808    0.432

Robust residual standard error: 1.281 
  (3 observations deleted due to missingness)
Multiple R-squared:  0.9923,	Adjusted R-squared:  0.9866 
Convergence in 1 IRWLS iterations

Robustness weights: 
 All 27 weights are ~= 1.
Algorithmic parameters: 
               bb        refine.tol           rel.tol         scale.tol 
        3.846e-04         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         3.704e-03         5.094e-12 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample     tuning.chi     tuning.psi         max.it       best.r.s 
           500             20             20             50              2 
      k.fast.s          k.max    maxit.scale      trace.lev            mts 
             1            200            200              0           1000 
    compute.rd fast.s.large.n 
             0           2000 
                  psi           subsampling                   cov 
            "optimal"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
> (sr2 <- summary(rm2))

Call:
lmrob(formula = y ~ x1 * x2 + x3 + x4 + x5, data = data2, offset = os, control = ctrl)
 \--> method = "MM"
Residuals:
     Min       1Q   Median       3Q      Max 
-1.52261 -0.51773  0.06925  0.38640  1.61986 

Coefficients: (3 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.002961   0.742168  -0.004    0.997
x1b         -1.132857   0.860835  -1.316    0.200
x1c          0.492904   0.985445   0.500    0.621
x2B          0.017959   0.922044   0.019    0.985
x2C          0.858031   0.914632   0.938    0.357
x2D          0.208510   0.969033   0.215    0.831
x3          -0.021632   0.215885  -0.100    0.921
x4          -0.079147   0.246574  -0.321    0.751
x5                 NA         NA      NA       NA
x1b:x2B      1.040529   1.096329   0.949    0.351
x1c:x2B            NA         NA      NA       NA
x1b:x2C      0.736944   1.162571   0.634    0.532
x1c:x2C            NA         NA      NA       NA
x1b:x2D      1.099090   1.110001   0.990    0.331
x1c:x2D     -1.371953   1.290375  -1.063    0.297

Robust residual standard error: 0.9728 
  (4 observations deleted due to missingness)
Multiple R-squared:  0.9923,	Adjusted R-squared:  0.989 
Convergence in 1 IRWLS iterations

Robustness weights: 
 All 38 weights are ~= 1.
Algorithmic parameters: 
               bb        refine.tol           rel.tol         scale.tol 
        3.846e-04         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         2.632e-03         4.369e-12 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample     tuning.chi     tuning.psi         max.it       best.r.s 
           500             20             20             50              2 
      k.fast.s          k.max    maxit.scale      trace.lev            mts 
             1            200            200              0           1000 
    compute.rd fast.s.large.n 
             0           2000 
                  psi           subsampling                   cov 
            "optimal"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
> 
> ## test Estimates, Std. Errors, ...
> nc <- names(coef(cm1))
> nc. <- setdiff(nc, "x5") # those who are "valid"
> stopifnot(exprs = {
+     all.equal(coef(cm0.),coef(cm0)[nc.])
+     all.equal(coef(cm1.),coef(cm1)[nc.])
+     all.equal(coef(cm2.),coef(cm2)[nc.])
+     all.equal(coef(cm1), coef(cm2))
+     all.equal(coef(rm1), coef(rm2))
+     all.equal(coef(sc0), coef(sr0))
+     all.equal(coef(sc1), coef(sr1))
+     all.equal(coef(sc2), coef(sr2))
+ })
> 
> ## test class "lm" methods that do not depend on weights
> meths1 <- c("family",
+             "formula",
+             "labels",
+             "model.matrix",
+             "na.action",
+             "terms")
> for (meth in meths1)
+     stopifnot(all.equal(do.call(meth, list(rm0)),
+                         do.call(meth, list(rm1))))
> 
> ## class "lm" methods that depend on weights
> ##                                      FIXME:
> meths2 <- c(#"AIC",
+             "alias",
+             #"BIC",
+             "case.names",
+             "coef",
+             "confint",
+             #"cooks.distance",
+             #"deviance",
+             "df.residual",
+             #"dfbeta",
+             #"dfbetas",
+             #"drop1",
+             "dummy.coef",
+             #"effects",
+             #"extractAIC",
+             #"hatvalues",
+             #"influence",
+             "kappa",
+             #"logLik",
+             #"model.frame", ## disable because of zero.weights attribute
+             "nobs",
+             "predict",
+                                         #"proj",
+                                         #"rstandard",
+                                         #"rstudent",
+                                         #"simulate",
+             ##"summary", ## see above
+             "variable.names",
+             ##"vcov",    ## see below
+             "weights")
> op <- options(warn = 1)# print immediately
> for (meth in meths2) {
+     cat(meth,":")
+     .SW. <- if(meth == "weights") suppressWarnings else identity # for suppressing
+     ## No weights defined for this object. Use type="robustness" ....
+     stopifnot(all.equal(do.call(meth, list(cm1)),
+                         do.call(meth, list(rm1))),
+               all.equal(do.call(meth, list(cm2)),
+ 		   .SW.(do.call(meth, list(rm2)))))
+     cat("\n")
+ }
alias :
case.names :
coef :
confint :
df.residual :
dummy.coef :
kappa :
nobs :
predict :
variable.names :
weights :
> options(op)# reverting
> 
> ## further tests:
> anova(rm1, update(rm1, ~ . - x4 - x5))
Robust Wald Test Table

Model 1: y ~ x1 * x2 + x3 + x4 + x5 + offset(os)
Model 2: y ~ x1 + x2 + x3 + x1:x2 + offset(os)
Largest model fitted by lmrob(), i.e. SM

  pseudoDf Test.Stat Df Pr(>chisq)
1       18                        
2       22  0.059442  1     0.8074
> anova(rm2, update(rm2, ~ . - x4 - x5))
Robust Wald Test Table

Model 1: y ~ x1 * x2 + x3 + x4 + x5
Model 2: y ~ x1 + x2 + x3 + x1:x2
Largest model fitted by lmrob(), i.e. SM

  pseudoDf Test.Stat Df Pr(>chisq)
1       23                        
2       27   0.10303  1     0.7482
> 
> stopifnot(all.equal(fitted(cm0), fitted(rm0)),
+           all.equal(fitted(cm1), fitted(rm1)),
+           all.equal(fitted(cm2), fitted(rm2)))
> 
> nd <- expand.grid(x1=letters[1:3], x2=LETTERS[1:4])
> set.seed(3)
> nd$x3 <- rnorm(nrow(nd))
> nd$x4 <- rnorm(nrow(nd))
> nd$x5 <- rnorm(nrow(nd)) # (*not* the sum x3+x4 !)
> nd$os <- nrow(nd):1
> wts   <- runif(nrow(nd))
> stopifnot(exprs = {
+     all.equal(predict(cm0, nd, interval="prediction"),
+               predict(rm0, nd, interval="prediction"))
+     all.equal(predict(cm1, nd, interval="prediction"),
+               predict(rm1, nd, interval="prediction"))
+     all.equal(predict(cm2, nd, interval="prediction"),
+               predict(rm2, nd, interval="prediction"))
+     all.equal(predict(cm0, nd, interval="prediction", weights=wts),
+               predict(rm0, nd, interval="prediction", weights=wts))
+     all.equal(predict(cm1, nd, interval="prediction", weights=wts),
+               predict(rm1, nd, interval="prediction", weights=wts))
+     all.equal(predict(cm2, nd, interval="prediction", weights=wts),
+               predict(rm2, nd, interval="prediction", weights=wts),
+               tolerance=1e-7)
+ })
There were 14 warnings (use warnings() to see them)
> 
> ## Padding can lead to differing values here
> ## so test only full rank part
> qrEQ <- function(m1, m2) {
+     q1 <- qr(m1)
+     q2 <- qr(m2)
+     r <- 1:q1$rank
+     stopifnot(q1$rank == q2$rank,
+               all.equal(q1$pivot, q2$pivot),
+               all.equal(q1$qraux[r],q2$qraux[r]),
+               all.equal(q1$qr[r,r], q2$qr[r,r]))
+ }
> qrEQ(cm0, rm0)
> qrEQ(cm1, rm1)
> qrEQ(cm2, rm2)
> 
> stopifnot(all.equal(residuals(cm0),             residuals(rm0)),
+           all.equal(residuals(cm1),             residuals(rm1)),
+           all.equal(residuals(cm2),             residuals(rm2)),
+           all.equal(resid(cm0, type="pearson"), resid(rm0, type="pearson")),
+           all.equal(resid(cm1, type="pearson"), resid(rm1, type="pearson")),
+           all.equal(resid(cm2, type="pearson"), resid(rm2, type="pearson")))
> 
> ## R 3.5.0:  vcov(*,  complete=TRUE)  new default ==> same NA's as coef()
> if(interactive()) withAutoprint({
+     op <- options(width = 130, digits = 2) # -->  vcov() rows fit on 1 line
+     vcov(cm0) # 'x5' is NA
+     vcov(cm2) # 'x5', 'x1c:2B', 'x1c:2C'  rows & columns are NA
+     options(op)
+ })
> 
> (no.C <- is.na(match("complete", names(formals(stats:::vcov.lm))))) ## temporary _FIXME_
[1] FALSE
> vcovC <- if(no.C) function(M, ...) vcov(M, complete=FALSE, ...) else vcov # (complete=TRUE)
> stopifnot(all.equal(vcov(cm0), vcovC(rm0), check.attributes=FALSE),
+           all.equal(vcov(cm1), vcovC(rm1), check.attributes=FALSE),
+           all.equal(vcov(cm2), vcovC(rm2), check.attributes=FALSE))
> 
> ## "clean":
> cln <- function(vc) structure(vc, weights=NULL, eigen=NULL)
> ## .vcov.avar1() is not recommended here, but also should work with singular / NA coef case:
> ok0 <- !is.na(coef(rm0))
> tools::assertWarning(verbose = TRUE, # on non-M1mac, there is a 2nd warning (not shown here!):
+ vr0.NA<- vcov(rm0, cov=".vcov.avar1", complete=NA))
Asserted warning: X'WX is almost singular. Consider using cov = ".vcov.w"
> vr0.T <- vcov(rm0, cov=".vcov.avar1", complete=TRUE)
> vr0.F <- vcov(rm0, cov=".vcov.avar1", complete=FALSE)
> stopifnot(identical(dim(vr0.NA), dim(vr0.T)),
+           identical(dim(vr0.F),  dim(vr0.T) - 1L), dim(vr0.F) == 14,
+           all.equal(cln(vr0.F), vr0.T[ok0,ok0], tol = 1e-15))
> 
> if(!no.C) {
+ vc0.T <- vcov(cm0, complete=TRUE)
+ vc0.F <- vcov(cm0, complete=FALSE)
+ }
> 
> ok1 <- !is.na(coef(rm1))
> ## cannot work because init/fit residuals are not of full length
> tools::assertError(vr1.NA<- vcov(rm1, cov=".vcov.avar1", complete=NA))
> tools::assertError(vr1.T <- vcov(rm1, cov=".vcov.avar1", complete=TRUE ))
> tools::assertError(vr1.F <- vcov(rm1, cov=".vcov.avar1", complete=FALSE))
> ## instead, must refit
> rm1. <- update(rm1, control = within(ctrl,  cov <- ".vcov.avar1"))
> vr1.NA<- vcov(rm1., complete=NA)
> vr1.T <- vcov(rm1., complete=TRUE)
> vr1.F <- vcov(rm1., complete=FALSE)
> 
> stopifnot(identical(vr1.F, vr1.NA), # in this case
+           identical(dim(vr1.F),  dim(vr1.T) - 3L), dim(vr1.F) == 12, isSymmetric(vr1.T),
+           identical(rownames(vr1.F),  rownames(vr1.T)[ok1]),
+           all.equal(cln(vr1.F), vr1.T[ok1,ok1], tol=1e-15))
> 
> if(FALSE) ## ERROR "exact singular" (probably *NOT* to fix, as TRUE/FALSE do work !)
+ vr2.NA<- vcov(rm2, cov=".vcov.avar1", complete=NA) # "almost singular" warning
> vr2.T <- vcov(rm2, cov=".vcov.avar1", complete=TRUE)
> vr2.F <- vcov(rm2, cov=".vcov.avar1", complete=FALSE)
> stopifnot(TRUE, # identical(dim(vr2.NA), dim(vr2.T)),
+ 	  identical(dim(vr2.F),  dim(vr2.T) - 3L), dim(vr2.F) == 12,
+ 	  identical(rownames(vr2.F), rownames(vr1.F)),
+ 	  identical(rownames(vr2.T), rownames(vr1.T)),
+ 	  all.equal(cln(vr2.F), vr2.T[ok1,ok1], tol=1e-15))
> 
> ## Hmm, the supposedly heteroscedastic-robust estimates *are* very different:
> all.equal(vcov(cm0), vcovC(rm0, cov = ".vcov.avar1"), check.attributes=FALSE) # rel.diff. 0.5367564
[1] "Mean relative difference: 0.5367564"
> if(FALSE) # does not make sense
+ all.equal(vcov(cm1), vcovC(rm1, cov = ".vcov.avar1"), check.attributes=FALSE)
> all.equal(vcov(cm2), vcovC(rm2, cov = ".vcov.avar1"), check.attributes=FALSE) # rel.diff. 0.5757642
[1] "Mean relative difference: 0.5757642"
> 
> 
> ## Null fits (rank(X)==0) are tested in NAcoef.R
> 
> ## testing weight=0 bug
> lmrob(y ~ x3, data, weights=weights)

Call:
lmrob(formula = y ~ x3, data = data, weights = weights)
 \--> method = "MM"
Coefficients:
(Intercept)           x3  
    18.7474       0.1751  

> 
> proc.time()
   user  system elapsed 
  0.315   0.082   0.450