tests/lmrob-data.R

In robustbase: Basic Robust Statistics

### lmrob()  with "real data" -----------------------

## testing functions:
source(system.file("xtraR/test-tools.R", package = "robustbase"))
(doExtras <- robustbase:::doExtras())
showProc.time()

library(robustbase)

##' short form to get "pure" robustness weights
rw <- function(fm) unname(weights(fm, type="robustness"))

set.seed(0)
data(salinity)
summary(m0.sali  <- lmrob(Y ~ . , data = salinity))
(A1 <- anova(m0.sali, Y ~ X1 + X3))
## -> X2 is not needed
(m1.sali  <- lmrob(Y ~ X1 + X3, data = salinity))
(A2 <- anova(m0.sali, m1.sali)) # the same as before
stopifnot(all.equal(A1[2,"Pr(>chisq)"],
		    A2[2,"Pr(>chisq)"], tolerance=1e-14))
anova(m0.sali, m1.sali, test = "Deviance")
## whereas 'X3' is highly significant:
m2 <- update(m0.sali, ~ . -X3)
(A3 <- anova(m0.sali, m2))
(A4 <- anova(m0.sali, m2, test = "Deviance"))
cX3 <- c(Estimate = -0.627327396, `Std. Error` = 0.15844971,
         `t value` = -3.9591577, `Pr(>|t|)` = 0.000584156)
stopifnot(all.equal(cX3, coef(summary(m0.sali))["X3",], tolerance = 1e-6))
showProc.time()


##  example(lmrob)
set.seed(7)
data(coleman)
summary( m1 <- lmrob(Y ~ ., data=coleman) )
stopifnot(identical(which(rw(m1) < 0.2), c(3L, 18L)))

if(FALSE) # to find out *why setting = "KS201x" fails
trace(lmrob.S, exit = quote({cat("coef:\n"); print(b$coefficients)}))
if(FALSE) # to find out via  setting = "KS201x" fails here in the *initial* estimate
debug(lmrob.S)

data(starsCYG)
  lmST <-    lm(log.light ~ log.Te, data = starsCYG)
(RlmST <- lmrob(log.light ~ log.Te, data = starsCYG, control=lmrob.control(trace = 1)))
summary(RlmST)
## Least Sq. w/ negative slope, where robust has slope ~= 2.2 :
stopifnot(exprs =  {
    coef(lmST)[["log.Te"]] < 0
    all.equal(coef(RlmST), c("(Intercept)" = -4.969, log.Te=2.253), tol = 1e-3)
    identical(which(rw(RlmST) < 0.01), as.integer( c(11,20,30,34) ))
})
showProc.time()
## ==> Now see that  "KS2011" and "KS2014" both break down -- and it is the fault of "lqq" *only* :
(RlmST.11 <- update(RlmST, control = lmrob.control("KS2011",                             trace= 1)))
(RlmST.14 <- update(RlmST, control = lmrob.control("KS2014",                             trace= 1)))
(RlmSTM11 <- update(RlmST, control = lmrob.control("KS2011", method="MM",                trace= 1)))
(RlmSTM14 <- update(RlmST, control = lmrob.control("KS2014", method="MM",                trace= 1)))
## using "biweight" instead of "lqq"  fixes the problem :
(RlmSTbM11 <- update(RlmST,control = lmrob.control("KS2011", method="MM", psi="biweight",trace= 1)))















(RlmSTbM14 <- update(RlmST,control = lmrob.control("KS2014", method="MM", psi="biweight",trace= 1)))
(RlmSTb.11 <- update(RlmST,control = lmrob.control("KS2011",              psi="biweight",trace= 1)))
(RlmSTb.14 <- update(RlmST,control = lmrob.control("KS2014",              psi="biweight",trace= 1)))
## NB: RlmST has component 'init.S' the others have "init" -- documented in ?lmrob.fit == ../man/lmrob.fit.Rd
R.ini.cf <- t(sapply(mget(ls(patt = "^RlmST")), function(OB) OB$init$coef))
R..cf    <- t(sapply(mget(ls(patt = "^RlmST")), coef))
cbind(R.ini.cf, R..cf) ##---> "lqq" is *NOT* robust enough here -- but "biweight" is !!

showProc.time()

options(digits = 5)# less platform dependence
## Directly look at init.S():
x.s <- model.matrix(~ log.Te, data = starsCYG)
y.s <- model.response(model.frame(log.light ~ log.Te, data = starsCYG))
ini.df <- lmrob.S(x.s, y.s, control=lmrob.control())
ini.11 <- lmrob.S(x.s, y.s, control=lmrob.control("KS2011"))
ini.14 <- lmrob.S(x.s, y.s, control=lmrob.control("KS2014"))
## but these are fine !! :
rbind(deflt = ini.df$coef,
      KS.11 = ini.11$coef,
      KS.14 = ini.14$coef)
##==> No, it is *not* the init.S()
ini.14$scale #  0.48144

## More clearly shows how M-estimate is converging to *WRONG* solution:
(RlmST.lqq <- update(RlmST, init=ini.14, control = lmrob.control(method="MM", psi="lqq", trace= 4)))
## --> break down
## The 10 largest residuals from the robust init. S-estim:
(i10 <- head(order(abs(residuals(ini.14)), decreasing=TRUE), 10))
residuals(ini.14)[i10]
## ==> and their weights for the different  psi() and their default (95% efficiency) tuning:
PSIs <- names(.Mchi.tuning.defaults)
sapply(PSIs, function(PSI)
    Mwgt(residuals(ini.14)[i10], cc = .Mpsi.tuning.defaults[[PSI]], psi=PSI))
## All MM:
RlmST.MM <- lapply(setNames(,PSIs), function(PSI)
    update(RlmST, init=ini.14, control = lmrob.control(method="MM", psi = PSI)))
cf.MM <- t(sapply(RlmST.MM, coef))
cf.MM[order(cf.MM[,1], cf.MM[,2]),] ## only 'bisquare' and 'optimal' are robust enough
showProc.time()


##=== Werner's analysis: Sensitivity curves for the most-left obs  =========================================
dd <- starsCYG
dd <- dd[order(dd[,"log.Te"]),] # ==> leverage points come first (and easier plotting)
(rr <- lmrob(log.light ~ log.Te, data = dd))
(rr14 <- update(rr, control = lmrob.control("KS2014")))

dd[1,2] # 6.05  will be replaced for sensitivity curve

leg.s <- c("default, biweight"
           ,"KS14, lqq"
           ,"KS14, biweight"
           ,"KS14, optimal"
           ,"KS14, Hampel"
           ,"KS14, GGW"
           ,"KS14, Welsh"
           )
nEst <- length(leg.s) # == number of estimators used below
nn <- length(y1 <- c(NA, seq(2,9, length=if(doExtras) 64 else 8)))
nCf <- length(coef(rr)) + 1 # +1: sigma
r.coef <- matrix(NA, length(y1), nEst*nCf)
t.d <- dd
oo <- options(warn = 1)
showProc.time()
## vary the left-most observation and fit all three
for (i in seq_along(y1)) {
  cat(sprintf("i=%3d, y1[i]=%11.6g: -- ", i, y1[i]))
  t.d[1,2] <- y1[i]
  ## the (old) default does not converge in 4 cases
  lr    <- update(rr,   data=t.d, control = lmrob.control(maxit=500))                ; cat(" 1")
  lr14  <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="lqq") )     ; cat(" 2")
  lr14b <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="biweight") ); cat(" 3")
  lr14o <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="optimal" ) ); cat(" 4")
  lr14h <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="hampel"  ) ); cat(" 5")
  lr14g <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="ggw"     ) ); cat(" 6")
  lr14a <- update(rr14, data=t.d, control = lmrob.control("KS2014", psi="welsh"   ) ); cat(" 7")
  r.coef[i,] <- c(coef(lr   ), sigma(lr),
                  coef(lr14 ), sigma(lr14),
                  coef(lr14b), sigma(lr14b),
                  coef(lr14o), sigma(lr14o),
                  coef(lr14h), sigma(lr14h),
                  coef(lr14g), sigma(lr14g),
                  coef(lr14a), sigma(lr14a))
  cat("\n")
}
showProc.time()
options(oo)

## cbind(y=y.1, r.coef)
## y1[1] = where the NA is
pMat <- function(j, main, x.legend, col = 1:8, lty=1:6, lwd = 2, ylab=NA, ...) {
    stopifnot(j %in% seq_len(ncol(r.coef)))
    matplot(y1, r.coef[, j], type="l", xlab = quote("varying obs." ~ ~ y[1]),
            ylab=ylab, main=main, col=col, lty=lty, lwd=lwd, ...)
    xx <- par("usr")[1:2]; yL <- .99* xx[1] + .01*xx[2]
    matpoints(yL, r.coef[1, j, drop=FALSE], pch = 20, col=col, lwd=lwd)
    abline(h =    r.coef[1, j, drop=FALSE], col = col, lwd=1, lty=3)
    legend(x.legend, leg.s, lty=lty, col=col, lwd=lwd, bty = "n")
    abline(v = dd[1,2], col=adjustcolor("tomato", 1/4)) # true y value
}

(jj0 <- nCf*(seq_len(nEst)-1L))
op <- {
    if(requireNamespace("sfsmisc", quietly=TRUE))
        sfsmisc::mult.fig(2)$old.par
    else
        par(mfrow = 2:1, mar = .1+ c(4,4,2,1), mgp = c(1.5, 0.6, 0))
}


pMat(j = 2+jj0, main = quote("slope" ~~ hat(beta[2])), "bottomleft")
pMat(j = 3+jj0, main = quote(hat(sigma)), "topleft")
par(op)
showProc.time()
## --------------------------------

set.seed(47)
data(hbk)
m.hbk <- lmrob(Y ~ ., data = hbk)
summary(m.hbk)
stopifnot(1:10 == which(rw(m.hbk) < 0.01))

data(heart)
summary(mhrt <- lmrob(clength ~ ., data = heart)) # -> warning 'maxit.scale=200' too small
stopifnot(8  == which(rw(mhrt) < 0.15),
          11 == which(0.61 < rw(mhrt) & rw(mhrt) < 0.62),
          c(1:7,9:10,12) == which(rw(mhrt) > 0.90))

iN <- c(3,5,7,11)
heartN <- heart; heartN[iN, "clength"] <- NA
lmN <- lm   (clength ~ ., data = heartN) # default na.action=na.omit
mhN <- lmrob(clength ~ ., data = heartN) # default na.action=na.omit
 # ==> everything just uses the n=8 complete obs
summary(mhN) # now *does* note the 4 omitted obs.
mhNex <- lmrob(clength ~ ., data = heartN, na.action=na.exclude)
summary(mhNex)
mhNx1 <- update(mhNex, ~ . - weight)
mhNx0 <- update(mhNex, ~ 1)
stopifnot(
    length(r.mNex <- resid(mhNex)) == nrow(heartN)
   ,
    iN == which(iNAr <- is.na(r.mNex))
   ,
    identical(iNAr, is.na(r.mN1 <- residuals(mhNx1)))
   ,
    identical(iNAr, is.na(r.mN0 <- residuals(mhNx0)))
)
showProc.time()

data(stackloss)
mSL <- lmrob(stack.loss ~ ., data = stackloss)
summary(mSL)

showProc.time()