R Under development (unstable) (2023-06-09 r84528 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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 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)) > vr0.NA<- vcov(rm0, cov=".vcov.avar1", complete=NA) # "almost singular" warning Warning messages: 1: In lf.cov(object, complete = complete, ...) : X'WX is almost singular. Consider using cov = ".vcov.w" 2: In lf.cov(object, complete = complete, ...) : .vcov.avar1: negative diag() fixed up; consider 'cov=".vcov.w."' instead > 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.37 0.12 0.48