-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
5 changed files
with
1,101 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,228 @@ | ||
| ##' Extract bread matrix from an lmrob fit | ||
| ##' | ||
| ##' This is part of a workaround for an issue in the robustbase code | ||
| ##' affecting sandwich covariance estimation. | ||
| ##' The issue in question is issue #6471, robustbase project on R-Forge. | ||
| ##' | ||
| ##' @title Bread method for lmrob objects | ||
| ##' @param x an lmrob object produced using an MM/SM estimator chain | ||
| ##' @param ... | ||
| ##' @return k by (k+1) matrix, with first column for scale estimate and rows, remaining cols for coefficients | ||
| ##' | ||
| ##' @author lrd author 2 | ||
| ##' @export | ||
| bread.lmrob <- function(x, ...) | ||
| { | ||
| stopifnot(is.list(ctrl <- x$control)) | ||
| if (!(!is.null(ctrl$method) && nchar(ctrl$method)<=2 && | ||
| substr(ctrl$method, nchar(ctrl$method),nchar(ctrl$method))=="M") ) | ||
| stop("bread.lmrob() supports only SM or MM estimates") | ||
|
|
||
| psi <- chi <- ctrl$psi | ||
| if (is.null(psi)) | ||
| stop("parameter psi is not defined") | ||
| stopifnot(is.numeric(c.chi <- ctrl$tuning.chi), is.numeric(c.psi <- ctrl$tuning.psi)) | ||
| r0 <- x$init$resid | ||
| r <- resid(x) | ||
| scale <- x$scale | ||
| xmat <- model.matrix(x) | ||
| bb <- 1/2 | ||
| n <- length(r) | ||
| stopifnot(n == length(r0), is.matrix(xmat), n == nrow(xmat)) | ||
| p <- ncol(xmat) | ||
| r.s <- r/scale | ||
| r0.s <- r0/scale | ||
| w <- robustbase::Mpsi(r.s, cc = c.psi, psi = psi, deriv = 1) | ||
| w0 <- robustbase::Mchi(r0.s, cc = c.chi, psi = chi, deriv = 1) | ||
| x.wx <- crossprod(xmat, xmat * w) | ||
| if (inherits(A <- tryCatch(solve(x.wx) * scale, error = function(e) e), | ||
| "error")) { | ||
| A <- tryCatch(solve(x.wx, tol = 0) * scale, error = function(e) e) | ||
| if (inherits(A, "error")) | ||
| { stop("X'WX is singular.") } else warning("X'WX is almost singular.") | ||
| } | ||
| ## At this point A has no sample size scaling, as in robustbase:::.vcov.avar1 | ||
| ## The lack of scaling there precisely compensates for the lack of scaling of the crossproduct | ||
| a <- A %*% (crossprod(xmat, w * r.s)/mean(w0 * r0.s)) | ||
| colnames(a) <- "sigma" | ||
| ## Now we restore sample size scaling to A | ||
| A <- n * A | ||
|
|
||
| cbind(a, A) | ||
| } | ||
| ##' | ||
| ##' | ||
| ##' Only SM or MM estimates supported | ||
| ##' | ||
| ##' @title Estfun method for lmrob objects | ||
| ##' @param x a fitted lmrob | ||
| ##' @param ... | ||
| ##' @return an estfun object, as in the sandwich package | ||
| ##' @author lrd author 2 | ||
| ##' @export | ||
| estfun.lmrob <- function(x, ...) | ||
| { | ||
| stopifnot(is.list(ctrl <- x$control)) | ||
| if (!(!is.null(ctrl$method) && nchar(ctrl$method)<=2 && | ||
| substr(ctrl$method, nchar(ctrl$method),nchar(ctrl$method))=="M") ) | ||
| stop("estfun.lmrob() supports only SM or MM estimates") | ||
|
|
||
| xmat <- model.matrix(x) | ||
| xmat <- naresid(x$na.action, xmat) | ||
| psi <- chi <- ctrl$psi | ||
| if (is.null(psi)) | ||
| stop("parameter psi is not defined") | ||
| stopifnot(is.numeric(c.chi <- ctrl$tuning.chi), is.numeric(c.psi <- ctrl$tuning.psi)) | ||
| r0 <- x$init$resid | ||
| r <- resid(x) | ||
| scale <- x$scale | ||
|
|
||
| n <- length(r) | ||
| stopifnot(n == length(r0), is.matrix(xmat), n == nrow(xmat)) | ||
| p <- ncol(xmat) | ||
| r0.s <- r0/scale | ||
| w0 <- robustbase::Mchi(r0.s, cc = c.chi, psi = chi) | ||
| Usigma <- scale(w0, center=TRUE, scale=FALSE) | ||
|
|
||
| r.s <- r/scale | ||
| w <- robustbase::Mpsi(r.s, cc = c.psi, psi = psi) | ||
|
|
||
| Ubeta <- w * xmat | ||
| rval <- cbind("sigma"=Usigma, Ubeta) | ||
| attr(rval, "assign") <- NULL | ||
| attr(rval, "contrasts") <- NULL | ||
| rval | ||
| } | ||
|
|
||
| ##' Overloading of sandwich::sandwich to accommodate non-square bread | ||
| ##' | ||
| ##' The sandwich package's sandwich function presumes the bread matrix | ||
| ##' to be symmetric. Obviously this won't do if the bread is rectangular but not | ||
| ##' square. | ||
| ##' | ||
| ##' This is part of a workaround for an issue in the robustbase code | ||
| ##' affecting sandwich covariance estimation. | ||
| ##' The issue in question is issue #6471, robustbase project on R-Forge. | ||
| ##' | ||
| ##' @title Sandwich estimate of covariance | ||
| ##' @param x a fitted model object, as in sandwich::sandwich | ||
| ##' @param bread. function or matrix, | ||
| ##' @param meat. function or matrix, as in sandwich::sandwich | ||
| ##' @param ... additional arguments to downstream methods, as in sandwich::sandwich | ||
| ##' @return matrix, bread %*% meat %*% t(bread) | ||
| ##' @author lrd author 2 | ||
| LRDsandwich <- function (x, bread. = sandwich::bread, meat. = sandwich::meat, ...) | ||
| { | ||
| if (is.list(x) && !is.null(x$na.action)) | ||
| class(x$na.action) <- "omit" | ||
| if (is.function(bread.)) | ||
| bread. <- bread.(x) | ||
| if (is.function(meat.)) | ||
| meat. <- meat.(x, ...) | ||
| n <- NROW(sandwich::estfun(x)) | ||
| ## the t() in the below is the only difference from sandwich::sandwich() | ||
| return(1/n * (bread. %*% meat. %*% t(bread.))) | ||
| } | ||
|
|
||
|
|
||
| .vcov.avar2 <- function(obj, x=obj$x, posdef.meth = c("posdefify","orig")) | ||
| { ## was .vcov.MM | ||
| stopifnot(is.list(ctrl <- obj$control)) | ||
| ## works only for MM & SM estimates: | ||
| ### if (!is.null(ctrl$method) && !ctrl$method %in% c('SM', 'MM')) | ||
| ### I replaced line above w/ 2 lines below for own reasons unrelated to bug fix -BH | ||
| if (!is.null(ctrl$method) && !(nchar(ctrl$method)==2 && | ||
| substr(ctrl$method, nchar(ctrl$method),nchar(ctrl$method))=="M") ) | ||
| stop('.vcov.avar2() supports only SM or MM estimates') | ||
| ## set psi and chi constants | ||
| psi <- chi <- ctrl$psi | ||
| if (is.null(psi)) stop('parameter psi is not defined') | ||
| stopifnot(is.numeric(c.chi <- ctrl$tuning.chi), | ||
| is.numeric(c.psi <- ctrl$tuning.psi)) | ||
|
|
||
| ## need (r0, r, scale, x, c.psi,c.chi, bb) | ||
| r0 <- obj$init$resid | ||
| r <- resid(obj) | ||
| scale <- obj$scale | ||
| if (is.null(x)) x <- model.matrix(obj) | ||
| bb <- 1/2 ## this is always 1/2 for S estimates by convention | ||
| ### --- start code from .vcov.MM --- | ||
| ## scaled residuals | ||
| n <- length(r) | ||
| stopifnot(n == length(r0), is.matrix(x), n == nrow(x)) | ||
| p <- ncol(x) | ||
| ## Next 2 lines added post-.vcov.MM, addressing #6471. | ||
| ## This assumes initial S-estimate solved sum( loss )/(n-p) == bb, | ||
| ## not mean( loss ) == bb as assumed in .vcov.avar1 | ||
| adj <- (n-p)/n | ||
| bb <- bb * adj | ||
| r.s <- r / scale # final scaled residuals | ||
| r0.s <- r0 / scale # initial scaled residuals | ||
| w <- Mpsi(r.s, cc = c.psi, psi = psi, deriv = 1) | ||
| w0 <- Mchi(r0.s, cc = c.chi, psi = chi, deriv = 1) | ||
| ## FIXME for multivariate y : | ||
| x.wx <- crossprod(x, x * w) | ||
| if(inherits(A <- tryCatch(solve(x.wx) * scale, | ||
| error=function(e)e), "error")) { | ||
| warning("X'WX is almost singular. Consider rather using cov = \".vcov.w\"") | ||
| A <- tryCatch(solve(x.wx, tol = 0) * scale, error=function(e)e) | ||
| if(inherits(A, "error")) | ||
| stop("X'WX is singular. Rather use cov = \".vcov.w\"") | ||
| } | ||
| a <- A %*% (crossprod(x, w * r.s) / mean(w0 * r0.s)) | ||
| w <- Mpsi( r.s, cc = c.psi, psi = psi) | ||
|
|
||
| ## 3) now the standard part (w, x, r0.s, n, A,a, c.chi, bb) | ||
| w0 <- Mchi(r0.s, cc = c.chi, psi = chi) | ||
| Xww <- crossprod(x, w*w0) | ||
| u1 <- A %*% crossprod(x, x * w^2) %*% (n * A) | ||
| u2 <- a %*% crossprod(Xww, A) | ||
| u3 <- A %*% tcrossprod(Xww, a) | ||
| u4 <- mean(w0^2 - bb^2) * tcrossprod(a) | ||
|
|
||
| ## list(cov = matrix((u1 - u2 - u3 + u4)/n, p, p), | ||
| ## wt = w / r.s, a = a) | ||
| ### --- end code from .vcov.MM --- | ||
| ret <- (u1 - u2 - u3 + u4)/n | ||
|
|
||
| ## this might not be a positive definite matrix | ||
| ## check eigenvalues (symmetric: ensure non-complex) | ||
| ev <- eigen(ret, symmetric = TRUE) | ||
| if (any(neg.ev <- ev$values < 0)) { ## there's a problem | ||
| posdef.meth <- match.arg(posdef.meth) | ||
| if(ctrl$trace.lev) | ||
| message("fixing ", sum(neg.ev), | ||
| " negative eigen([",p,"])values") | ||
| Q <- ev$vectors | ||
| switch(posdef.meth, | ||
| "orig" = { | ||
| ## remove negative eigenvalue: | ||
| ## transform covariance matrix into eigenbasis | ||
| levinv <- solve(Q) | ||
| cov.eb <- levinv %*% ret %*% Q | ||
| ## set vectors corresponding to negative ev to zero | ||
| cov.eb[, neg.ev] <- 0 | ||
| ## cov.eb[cov.eb < 1e-16] <- 0 | ||
| ## and transform back | ||
| ret <- Q %*% cov.eb %*% levinv | ||
| }, | ||
| "posdefify" = { | ||
| ## Instead of using require("sfsmisc") and | ||
| ## ret <- posdefify(ret, "someEVadd",eigen.m = ev,eps.ev = 0) | ||
| lam <- ev$values | ||
| lam[neg.ev] <- 0 | ||
| o.diag <- diag(ret)# original one - for rescaling | ||
| ret <- Q %*% (lam * t(Q)) ## == Q %*% diag(lam) %*% t(Q) | ||
| ## rescale to the original diagonal values | ||
| ## D <- sqrt(o.diag/diag(ret)) | ||
| ## where they are >= 0 : | ||
| D <- sqrt(pmax.int(0, o.diag)/diag(ret)) | ||
| ret[] <- D * ret * rep(D, each = p) ## == diag(D) %*% m %*% diag(D) | ||
| }, | ||
| stop("invalid 'posdef.meth': ", posdef.meth)) | ||
| } | ||
| attr(ret,"weights") <- w / r.s | ||
| attr(ret,"eigen") <- ev | ||
| ret | ||
| }## end{.vcov.avar2} | ||
|
|
Oops, something went wrong.