R/linreg.R

In Sidryd/lab4sidjac: Linear regression

#' Fitting linear regression
#' @description \code{linreg} is used to fit linear regression models.
#' @param formula an object of class \code{\link[stats]{formula}}.
#' @param data a data frame containing all the variables in the regression model stated in the argument \code{formula}.
#' @return \code{linreg} return a named list of class "linreg".
#' \itemize{
#' \item{\code{call_arg}}{ return the arguments used to perform the linear regression}
#' \item{\code{y}}{ a vector with the response variable data used}
#' \item{\code{X}}{ the dependent variables expressed with \code{\link[stats]{model.matrix}}}
#' \item{\code{beta_hat}}{ a one column matrix with the coefficients}
#' \item{\code{sigma_2_hat}}{ the residual variance in a 1x1 matrix}
#' \item{\code{y_hat}}{ a one column matrix with the fitted values}
#' \item{\code{e_hat}}{ a one column matrix with the residuals}
#' \item{\code{var_beta_hat}}{ the covariance matrix for the coefficients}
#' \item{\code{t_beta}}{ a matrix with t-values for each coefficient}
#' \item{\code{p_value}}{ a matrix with p-values for each coefficient}
#' \item{\code{data}}{ return the data.frame stated in the argument \code{data}}}
#' @examples
#' linreg(formula = Petal.Length ~ Sepal.Length + Sepal.Width, data = iris)
#' @references https://en.wikipedia.org/wiki/Linear_regression
#' @export

# Main function
linreg <- function(formula, data){
  call_arg <- match.call()
  y <- data[,all.vars(formula)[1]]
  X <- stats::model.matrix(object = stats::formula(formula)[c(1,3)], data = data)
  beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
  y_hat <- X %*% beta_hat
  e_hat <- y - y_hat
  df <- nrow(data) - ncol(X)
  sigma_2_hat <- (t(e_hat) %*% e_hat) / df
  var_beta_hat <- sigma_2_hat[1] * solve(t(X) %*% X)
  t_beta <- beta_hat / sqrt(diag(var_beta_hat))
  p_value <- 2*stats::pt(abs(t_beta), 147, lower.tail = FALSE)
  res <- list(call_arg = call_arg, y = y, X = X, beta_hat = beta_hat,
              y_hat = y_hat, e_hat = e_hat, df = df, sigma_2_hat = sigma_2_hat,
              var_beta_hat = var_beta_hat, t_beta = t_beta, p_value = p_value,
              data = data)
  class(res) <- "linreg"
  invisible(res)
}