R/theilsen.R

In deming: Deming, Theil-Sen, Passing-Bablock and Total Least Squares Regression

# Function for Theil-Sen regression
theilsen <- function(formula, data, subset, weights, na.action, conf=.95,
                  nboot=0, symmetric=FALSE, 
                  eps=sqrt(.Machine$double.eps),
                  x=FALSE, y=FALSE, model=TRUE) {
    Call <- match.call()

    # create a call to model.frame() that contains the formula (required)
    #  and any other of the relevant optional arguments
    # then evaluate it in the proper frame
    indx <- match(c("formula", "data", "weights", "subset", "na.action"),
                  names(Call), nomatch=0) 
    if (indx[1] ==0) stop("A formula argument is required")
    temp <- Call[c(1,indx)]  # only keep the arguments we wanted
    temp[[1]] <- as.name('model.frame')  # change the function called
    mf <- eval(temp, parent.frame())
    Terms <- terms(mf)

    X <- model.matrix(Terms, mf)
    if (!attr(Terms, "intercept")) stop ("an intercept is required")
    if (ncol(X) != 2) 
        stop("Theil-Sen regression requires a single predictor variable")

    Y <- model.response(mf, type="numeric")
    n <- length(Y)
    if (is.null(Y))
        stop ("a response variable is required")
    if (n <3) stop("less than 3 non-missing observations in the data")

    wt <- model.weights(mf)
    if (is.null(wt)) wt <- rep(1.0, n)
    else {
        if (conf >0 && nboot==0) 
            stop("confidence intervals for weighted data require a bootstrap")
    }

    if (!is.numeric(eps) || eps <=0) stop("invalid value for eps")
    if (!is.logical(symmetric)) stop("invalid value for symmetric argument")
    if (conf <0 || conf >=1) stop("invalid confidence interval limit")
    
    fit <- theilsen.fit(X[,2], Y, wt, symmetric, conf*(nboot==0), eps)
    names(fit$coefficients) <-  dimnames(X)[[2]]
    yhat <- fit$coefficients[1] + X[,2]*fit$coefficients[2]
    fit$residuals <- Y-yhat

    if (nboot > 0) {
        # Set up for the wild bootstrap, which add random noise to
        #  the original x,y values
        # "residuals" are orthagonal to the fitted line
        # cos(theta) = 1/sqrt(1+ slope^2),  sin(theta) = slope/sqrt(1+slope^2)
        fcoef <- fit$coefficients
        d <- sqrt(1+ fcoef[2]^2)
        u <- (X[2] + fcoef[2]*(Y - fcoef[1]))/d  #projection onto the fit
        px <- u/d; py <- fcoef[1] + fcoef[2]*u/d
        resid <- list(x= X[2]- px, y= Y- py)
 
        wild <- function(data, resid) {
            n <- length(data$x)
            rb <- matrix(rbinom(2*n, 1, prob=(sqrt(5)+1)/sqrt(20)), ncol=2)
            newx <- data$x + resid$x*ifelse(rb[,1]==1, -(sqrt(5)-1)/2,
                                             (sqrt(5)+1)/2)
            newy <- data$y + resid$y*ifelse(rb[,2]==1, -(sqrt(5)-1)/2,
                                             (sqrt(5)+1)/2)
            list(x=newx, y=newy, w=data$w)
        }
        dfun <- function(data) {  #will be called by boot function
            theilsen.fit(data$x, data$y, data$w, symmetric, conf=0, 
                           eps)$coefficients
         }
        boot.out <- boot(list(x=X[,2], y=Y, w= wt), dfun, nboot, 
                         sim="parametric",
                         ran.gen=wild, mle=resid)
        temp1 <- boot.ci(boot.out, type="perc", index=1, conf=conf)
        temp2 <- boot.ci(boot.out, type="perc", index=2, conf=conf)
        fit$boot <-boot.out
        fit$ci   <- rbind(temp1$perc[,4:5], temp2$perc[,4:5])
        fit$var <- var(boot.out$t)
    }

 #   dimnames(fit$ci) <- list(names(fit$coef), 
#                             paste(c("lower", "upper"), format(conf)))
    fit$n <- n
    fit$conf <- conf
    if (x) fit$x <- X
    if (y) fit$y <- Y
    if (model) fit$model <- mf
    
    na.action <- attr(mf, "na.action")
    if (length(na.action)) fit$na.action <- na.action
    fit$terms <- Terms
    class(fit) <- c("theilsen") 
    fit$call <- Call
    fit
}


print.theilsen <- function(x, ...) {
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    cat("n=", x$n)
    if (length(x$na.action))
        cat("  (", naprint(x$na.action), ")\n", sep='')
    else cat("\n")

    temp <- matrix(x$coefficients, 2,1)
    tname <- "Coefficient"
    if (!is.null(x$var)) {
        temp <- cbind(temp, sqrt(diag(x$var)))
        tname <- c(tname, "Std err")
    }
    if (!is.null(x$ci)){
        temp <- cbind(temp, x$ci)
        tname <- c(tname, paste(c("lower", "upper"), format(signif(x$conf,2))))
    }
    dimnames(temp) <- list(names(x$coefficients), tname)

    print(temp, ...)
    invisible(x)
    }