R/slope.R

In rtrend: Trend Estimating Tools

#' slope
#'
#' * `slope`     : linear regression slope
#' * `slope_p`   : linear regression slope and p-value
#' * `slope_mk`  : mann kendall Sen's slope and p-value
#' * `slope_sen` : same as `slope_mk`, but with no p-value
#' * `slope_boot`: bootstrap slope and p-value
#'
#' @param y vector of observations of length n, or a matrix with n rows.
#' @param x vector of predictor of length n, or a matrix with n rows.
#' @param fast Boolean. If true, [stats::.lm.fit()] will be used, which is 10x
#' faster than [stats::lm()].
#' @param ... ignored.
#' 
#' @return
#' - `slope`  : linear regression coefficient
#' - `pvalue` : `p-value <= 0.05`` means that corresponding `slope` is significant.
#' - `sd`     : `Std. Error`
#'
#' For `slope_boot`, slope is estimated in many times. The lower, mean, upper
#' and standard deviation (sd) are returned.
#'
#' @examples
#' y <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
#' r    <- slope(y)
#' r_p  <- slope_p(y)
#' r_mk <- slope_mk(y)
#' r_boot <- slope_boot(y)
#' @importFrom stats .lm.fit pt
#' @export
slope <- function(y, x, ...){
    # TODO: add tests for slop
    if (!is.matrix(y)) y <- as.matrix(y)
    n <- nrow(y)

    if (missing(x)) x <- as.matrix(1:n)
    if (!is.matrix(x)) x <- as.matrix(x)

    I_bad <- which(!is.finite(y)) # NA or Inf

    if (length(I_bad) > 0) {
        y <- y[-I_bad, , drop = FALSE]
        x <- x[-I_bad, , drop = FALSE]
    }
    if (length(y) <= 1) return(c(slope = NA_real_))

    slope = qr.solve(cbind(x*0+1, x) , y)[2, ]
    if (length(slope) == 1) names(slope) = "slope"
    slope
}

#' @rdname slope
#' @export
slope_p <- function(y, x, fast = TRUE){
    if (!is.matrix(y)) y <- as.matrix(y)
    n <- nrow(y)

    if (missing(x)) x <- as.matrix(1:n)
    if (!is.matrix(x)) x <- as.matrix(x)

    I_bad <- which(!is.finite(y)) # NA or Inf

    if (length(I_bad) > 0) {
        y <- y[-I_bad, , drop = FALSE]
        x <- x[-I_bad, , drop = FALSE]
    }
    if (length(y) <= 1) return(c(slope = NA_real_, pvalue = NA_real_))

    # pvalue: the smaller, the better
    if (fast) {
        l <- .lm.fit(x = cbind(1, x), y = y)
        coefficients <- summary_lm(l)
    } else {
        l <- lm(y ~ x)
        coefficients <- summary(l)$coefficients
    }
    coefficients[2, c(1, 4, 2)] %>% set_names(c("slope", "pvalue", "sd"))
}

#' @rdname slope
#' @export
slope_sen_r <- function(y, x = seq_along(y), ...) {
    n = length(x)
    V <- rep(NA, times = (n^2 - n) / 2)
    k = 0
    for (i in 2:n) {
        for (j in 1:(i - 1)) {
            # for (j in 1:(n - 1)) {
            k = k + 1
            V[k] = (y[i] - y[j]) / (x[i] - x[j])
        }
    }
    median(na.omit(V))
}

#' @rdname slope
#' @export
slope_mk <- function(y, x = NULL, ...){
    mkTrend(y, x)[c("slp", "pval")] %>% set_names(c("slope", "pvalue"))
}

#' @param slope_FUN one of [slope()], [slope_p()], [slope_mk()]
#' @param times The number of bootstrap replicates.
#' @param alpha significant level, defalt 0.1
#'
#' @inheritParams boot::boot
#' @inheritParams base::set.seed
#'
#' @rdname slope
#' @importFrom boot boot
#' @importFrom matrixStats colQuantiles colSds
#' @export
slope_boot <- function(y, x = NULL, slope_FUN = slope, times = 100, alpha = 0.1, seed, ...) {
    if (!missing(seed)) set.seed(seed)

    x <- x %||% seq_along(y)
    FUN <- function(y0, indices) {
        y2 <- y[indices]
        x2 <- x[indices]
        slope_FUN(y2, x2)
    }
    b <- boot(y, FUN, R = times)

    probs  <- c(alpha/2, 0.5, 1 - alpha/2)
    b$coef <- colQuantiles(as.matrix(b$t), probs = probs, drop = FALSE) %>%
        set_rownames(names(b$t0)) %>%
        set_colnames(c("lower", "mean", "upper")) %>%
        cbind(sd = colSds(b$t))
    b$coef
}