# R/conRLM_fit.R In restriktor: Restricted Statistical Estimation and Inference for Linear Models

```conRLM_fit <- function (x, y, weights, ..., w = rep(1, nrow(x)), init = "ls",
psi = psi.bisquare, scale.est = c("MAD", "Huber", "proposal 2"),
k2 = 1.345, method = c("M", "MM"),
wt.method = c("inv.var", "case"),
maxit = 20, acc = 1e-04, test.vec = "resid",
lqs.control = NULL, Amat, bvec, meq,
tol = sqrt(.Machine\$double.eps)) {

irls.delta <- function(old, new) sqrt(sum((old - new)^2)/max(1e-20,
sum(old^2)))
irls.rrxwr <- function(x, w, r) {
w <- sqrt(w)
max(abs((matrix(r * w, 1L, length(r)) %*% x)/sqrt(matrix(w,
1L, length(r)) %*% (x^2))))/sqrt(sum(w * r^2))
}
o <- sort.list(abs(x))
x <- abs(x)[o]
w <- w[o]
p <- cumsum(w)/sum(w)
n <- sum(p < 0.5)
if (p[n + 1L] > 0.5)
x[n + 1L]/0.6745
else (x[n + 1L] + x[n + 2L])/(2 * 0.6745)
}
method <- match.arg(method)
wt.method <- match.arg(wt.method)
nmx <- deparse(substitute(x))
if (is.null(dim(x))) {
x <- as.matrix(x)
colnames(x) <- nmx
}
else x <- as.matrix(x)
if (is.null(colnames(x)))
colnames(x) <- paste("X", seq(ncol(x)), sep = "")
if (qr(x)\$rank < ncol(x))
stop("'x' is singular: singular fits are not implemented in 'rlm'")
if (!(any(test.vec == c("resid", "coef", "w", "NULL")) ||
is.null(test.vec)))
stop("invalid 'test.vec'")
xx <- x
yy <- y
if (!missing(weights)) {
if (length(weights) != nrow(x))
stop("length of 'weights' must equal number of observations")
if (any(weights < 0))
stop("negative 'weights' value")
if (wt.method == "inv.var") {
fac <- sqrt(weights)
y <- y * fac
x <- x * fac
wt <- NULL
}
else {
w <- w * weights
wt <- weights
}
}
else wt <- NULL
if (method == "M") {
scale.est <- match.arg(scale.est)
if (!is.function(psi))
psi <- get(psi, mode = "function")
arguments <- list(...)
if (length(arguments)) {
pm <- pmatch(names(arguments), names(formals(psi)),
nomatch = 0L)
if (any(pm == 0L))
warning("some of ... do not match")
pm <- names(arguments)[pm > 0L]
formals(psi)[pm] <- unlist(arguments[pm])
}
if (is.character(init)) {
temp <- if (init == "ls")
lm.wfit(x, y, w, method = "qr")
else if (init == "lts") {
if (is.null(lqs.control))
lqs.control <- list(nsamp = 200L)
do.call("lqs", c(list(x, y, intercept = FALSE),
lqs.control))
}
else stop("'init' method is unknown")
coef <- temp\$coefficient
resid <- temp\$residuals
}
else {
if (is.list(init))
coef <- init\$coef
else coef <- init
resid <- drop(y - x %*% coef)
}
}
else if (method == "MM") {
scale.est <- "MM"
psi <- psi.bisquare

if (any(bvec != 0)) {
stop("Restriktor ERROR: the rhs may only contain zeros (for now).")
}

if (meq > 0L) {
Dmat <- crossprod(x)
dvec <- t(x) %*% y
out.qp <- solve.QP(Dmat, dvec, t(Amat[1:meq, , drop = FALSE]),
bvec[1:meq], meq = meq)\$solution
out.qp[abs(out.qp) < tol] <- 0L
idx.qp <- !out.qp %in% 0L
temp <- do.call("lqs",
c(list(x = x[ , idx.qp, drop = FALSE], y, intercept = FALSE,
method = "S", k0 = 1.54764), lqs.control))
coef  <- temp\$coefficients
resid <- temp\$residuals
scale <- temp\$scale
} else {
temp <- do.call("lqs", c(list(x, y, intercept = FALSE,
method = "S", k0 = 1.548), lqs.control))
coef <- temp\$coefficients
resid <- temp\$residuals
}
if (length(arguments <- list(...)))
if (match("c", names(arguments), nomatch = 0L)) {
c0 <- arguments\$c
if (c0 > 1.54764)
formals(psi)\$c <- c0
else warning("'c' must be at least 1.548 and has been ignored")
}
scale <- temp\$scale
} else {
stop("'method' is unknown")
}

done <- FALSE
conv <- NULL
n1 <- (if (is.null(wt))
nrow(x)
else sum(wt)) - ncol(x)
theta <- 2 * pnorm(k2) - 1
gamma <- theta + k2^2 * (1 - theta) - 2 * k2 * dnorm(k2)
if (scale.est != "MM")
scale <- if (is.null(wt))
for (iiter in 1L:maxit) {
if (!is.null(test.vec))
testpv <- get(test.vec)
if (scale.est != "MM") {
scale <- if (scale.est == "MAD")
if (is.null(wt))
median(abs(resid))/0.6745
else if (is.null(wt))
sqrt(sum(pmin(resid^2, (k2 * scale)^2))/(n1 *
gamma))
else sqrt(sum(wt * pmin(resid^2, (k2 * scale)^2))/(n1 *
gamma))
if (scale == 0) {
done <- TRUE
break
}
}
w <- psi(resid/scale)
if (!is.null(wt))
w <- w * weights

#############################################################
W <- sqrt(w)
fit <- con_solver_rlm(X = x * W, y = y * W,
Amat = Amat, bvec = bvec,
meq = meq, maxit = maxit,
absval = tol)

coef <- fit\$solution
resid <- drop(y - x %*% coef)
iact <- fit\$iact
#############################################################

#temp <- lm.wfit(x, y, w, method = "qr")
#coef <- temp\$coefficients
#resid <- temp\$residuals
if (!is.null(test.vec))
convi <- irls.delta(testpv, get(test.vec))
else convi <- irls.rrxwr(x, w, resid)
conv <- c(conv, convi)
done <- (convi <= acc)
if (done)
break
}
if (!done)
warning(gettextf("'rlm' failed to converge in %d steps",
maxit), domain = NA)
fitted <- drop(xx %*% coef)
cl <- match.call()
cl[[1L]] <- as.name("rlm")
fit <- list(coefficients = coef, residuals = yy - fitted,
wresid = resid, effects = temp\$effects, rank = temp\$rank,
fitted.values = fitted, assign = temp\$assign, qr = temp\$qr,
df.residual = NA, w = w, s = scale, psi = psi, k2 = k2,
weights = if (!missing(weights)) weights, conv = conv,
converged = done, x = xx, call = cl,
Amat = Amat, bvec = bvec, meq = meq, iact = iact)
class(fit) <- c("conRLM", "rlm", "lm")
fit
}
```

## Try the restriktor package in your browser

Any scripts or data that you put into this service are public.

restriktor documentation built on Feb. 25, 2020, 5:08 p.m.