R/summary_lm.R

In rtrend: Trend Estimating Tools

#' summary_lm
#' 
#' summary method for class ".lm.fit".. It's 200 times faster than traditional \code{lm}.
#' 
#' @param obj Object returned by \code{\link{.lm.fit}}.
#' @param ... ignored
#' 
#' @return a p x 4 matrix with columns for the estimated coefficient, its standard error, 
#' t-statistic and corresponding (two-sided) p-value. Aliased coefficients are omitted.
#' 
#' @example R/examples/ex-summary_lm.R
#' @export
summary_lm <- function(obj, ...){
    z <- obj
    p <- z$rank

    Qr <- z$qr
    n <- nrow(Qr)
    rdf <- n - p#df.residual

    r <- z$residuals
    # f <- z$fitted.values
    # mss <- sum((f - mean(f))^2)
    rss <- sum(r^2)
    resvar <- rss/rdf

    p1 <- 1L:p
    R <- chol2inv(Qr[p1, p1, drop = FALSE])
    se <- sqrt(diag(R) * resvar)
    est <- z$coefficients[z$pivot[p1]]
    tval <- est/se
    # ans <- z[c("call", "terms", if (!is.null(z$weights)) "weights")]
    # ans$residuals <- r
    coefficients <- cbind(est, se, tval, 2 * pt(abs(tval),
                                                rdf, lower.tail = FALSE))
    dimnames(coefficients) <- list(names(z$coefficients)[z$pivot[p1]],
                                   c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
    return(coefficients)
    # df.int <- ifelse(attr(z$terms, "intercept"), 1L, 0L)
    # ans$r.squared <- mss/(mss + rss)
    # ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n - df.int)/rdf)
    # ans$fstatistic <- c(value = (mss/(p - df.int))/resvar, numdf = p - df.int, dendf = rdf)
    # ans
}