R/linreg.R

In eliscl/lab4: Estimates linear regression.

#'@title Function to estimate (multiple) linear regression.
#'@param formula an object of class "formula": a symbolic description of the model to be fitted.
#'@param data the dataframe containing the variables in the model.
#'@source \url{https://en.wikipedia.org/wiki/Linear_regression}
#'@export

linreg <- function(formula, data) {
  ## PREPARING DATA FOR ESTIMATION
  # Extracting a matrix with the independent variables.
  X <- model.matrix(formula, data = data)
  # Extracting a matrix with the dependent variable.
  y <- all.vars(formula)[1]
  y <- as.matrix(data[y])

  ## ESTIMATING PARAMETERS
  # Beta coefficients
  coefficients <- c("Coefficients" = solve(t(X) %*% X) %*% t(X) %*% y)
  # Predicted values
  fitted_values <- X %*% coefficients
  # Residuals
  residuals <- y - fitted_values
  # Degrees of freedom
  n <- nrow(data) # nr. of observations
  p <- ncol(X) # nr. of variables
  df <- n - p # degrees of freedom
  # var(betas)
  var <- 1 / df * as.numeric(t(residuals) %*% residuals) * solve(t(X) %*% X) # Variance-covariance matrix
  #var <- var(t(res)%*%res) / df
  var <- diag(var)
  # standard errors
  se <- sqrt(diag(var))
  se <- diag(se)
  # t-values
  t_value <- coefficients / sqrt(var)
  # p-value
  p_value <-2 * pt(abs(t_value), df, lower.tail = F)

  # OUTPUT
  # Merging all information
  # Adding the output to a S3 object.
  output <-
    structure(
      list(
        "Formula" = formula,
        "Coefficients" = coefficients,
        "Residuals" = as.vector(residuals),
        "Std.e." = se,
        "Fitted values" = fitted_values,
        "T-values" = t_value,
        "P-values" = p_value,
        "df" = df,
        "call" = match.call()
      ),
      class = "linreg"
    )
  names(output$Coefficients) <- colnames(X)
  return(output)
}