R/linmodEst.R

In janek01/testPackage:

#' Linear regression
#'
#' Runs an OLS regression not unlike \code{\link{lm}}
#'
#' @param y response vector (1 x n)
#' @param X covariate matrix (p x n) with no intercept
#'
#' @return A list with 4 elements: coefficients, vcov, sigma, df
#'
#' @examples
#' data(mtcars)
#' X <- as.matrix(mtcars[, c('cyl', 'disp', 'hp')])
#' y <- mtcars[, 'mpg']
#' linreg(y, X)
#'
#' @export
#'

linmodEst <- function(x, y) {
    ## CC: crossprod or a QR decomposition (as in the original version) are more efficient
    coef <- solve(t(x) %*% x) %*% t(x) %*% y
    print(coef)
    ## degrees of freedom and standard deviation of residuals
    df <- nrow(x) - ncol(x)
    sigma2 <- sum((y - x %*% coef)^2)/df
    ## compute sigma^2 * (x’x)^-1
    vcov <- sigma2 * solve(t(x) %*% x)
    colnames(vcov) <- rownames(vcov) <- colnames(x)
    list(coefficients = coef, vcov = vcov, sigma = sqrt(sigma2), df = df)
}

# Run: data(cats, package = 'MASS') linmodEst(cbind(1, cats$Bwt), cats$Hwt)

linmod <- function(x, ...)
  UseMethod("linmod")

#核心计算函数
linmod.default <- function(x, y, ...) {
  x <- as.matrix(x)
  y <- as.numeric(y)
  est <- linmodEst(x, y)
  est$fitted.values <- as.vector(x %*% est$coefficients)
  est$residuals <- y - est$fitted.values
  est$call <- match.call()
  class(est) <- "linmod"
  return(est)
}

#定义输出
print.linmod <- function(x, ...) {
  cat("Call:\n")
  print(x$call)
  cat("\nCoefficients:\n")
  print(x$coefficients)
}

#定义输出格式
linmod.formula <- function(formula, data = list(), ...) {
  mf <- model.frame(formula = formula, data = data)
  x <- model.matrix(attr(mf, "terms"), data = mf)
  y <- model.response(mf)
  est <- linmod.default(x, y, ...)
  est$call <- match.call()
  est$formula <- formula
  return(est)
}

#Run: linmod(Hwt ~ - 1 + Bwt * Sex, data = cats)


#TEST:
#devtools::use_testthat()