R/lm_fit.R

In leyaozh/lm.hw4: Linear Regression

#'lm_fit
#'
#'Fit a linear regression model.
#'@param x a n * p matrix of predictor variables.
#'@param y a length n vector(q = 1) or a n * q matrix of response variables.
#'@param add.intercept logical. If TRUE, add intercept term to the linear regression model.
#'@param method the method to fit the model(only method "qr" and "inv" is supported). Default is qr.
#'@param tol the tolerance for detecting linear dependencies in the columns of x. Default is 1e-7.
#'
#'@return
#' lm_fit returns a list containing following response values.
#'\item{coefficients}{q named vectors of coefficients.}
#'\item{residuals}{q vectors of residuals, which equals response minus fitted values.}
#'\item{rank}{the numeric rank of the fitted linear model.}
#'\item{fitted.values}{q vectors of the fitted mean values.}
#'\item{df.residual}{the degrees of freedom of residuals.}
#'\item{method}{the method to be used for fitting the model.}
#'
#' The following response values are returned and printed only when method = "qr".
#'\item{effects}{q vectors of orthogonal single-df effects.}
#'\item{qr}{the QR decomposition of matrix predictor variables.}
#' The functions lm_summary and lm_anova are used to obtain a summary and analysis of variance table of the results.
#'
#'@examples
#'n = 10; p = 5; q = 2;
#'x = matrix(rnorm(n * p), n, p) # no intercept
#'y1 = rnorm(n)
#'y2 = matrix(rnorm(n * q), n, q)
#'
#'# no intercept, using method "qr"
#'z = lm_fit(x = x, y = y1)
#'# no intercept, using method "inv"
#'z = lm_fit(x = x, y = y1, method ="inv")
#'# with intercept, using method "qr"
#'z = lm_fit(x = x, y = y2, add.intercept = TRUE)
#'# with intercept, using method "inv"
#'z = lm_fit(x = x, y = y2, add.intercept = TRUE, method = "inv")
#'
#'@export
#'

lm_fit <- function(x, y, add.intercept = FALSE, method = "qr", tol = 1e-07){
  nrx = nrow(x)
  ncx = ncol(x)
  y = as.matrix(y)
  nry = nrow(y)
  ncy = ncol(y)

  if (is.null(nrx)){
    stop("'x' must be a matrix")
  }
  if (!is.null(nry) && nrx != nry){
    stop("incompatible dimensions")
  }
  if(add.intercept == TRUE){
    x = cbind(rep(1, nrow(x)), x)
    ncx = ncx + 1
  }
  if (qr(x)$rank < ncx){
    stop("singular fit encountered")
  }
  if (method != "inv" && method != "qr"){
    stop("only support 'qr' and 'inv' method")
  }

  # if add.intercept == TRUE, add an intercept term to the model.
  if (add.intercept){
    name = c("(Intercept)",paste0("x", 1 : (ncx-1)))
  }else{
    name = paste0("x", 1 : ncx)
  }

  # estimate parameters using inverse matrix.
  if (method == "inv"){
    z <- NULL
    coefficients = solve(t(x) %*% x) %*% t(x) %*% y
    fitted.values = x %*% coefficients
    residuals = y - fitted.values
    df.residuals = nrx - ncx

    if (ncy == 1){
      coefficients = drop(coefficients)
      names(coefficients) = name
      fitted.values = drop(fitted.values)
      residuals = as.vector(residuals)
    }else{
      rownames(coefficients) = name
    }

    z$add.intercept = add.intercept
    z$coefficients = coefficients
    z$residuals = residuals
    z$rank = rank(x)
    z$fitted.values = fitted.values
    z$df.residual = df.residuals
    z$method = method
  }

  # estimate parameters using qr method.
  if (method == "qr"){
    z <- NULL
    qr = qr(x, tol) #QR decomposition of X, X = QR
    effects = qr.qty(qr(x,tol),y)
    coefficients = solve.qr(qr,y) # coefficients = R^(-1)Q^TY
    fitted.values = x %*% coefficients
    residuals = y - fitted.values
    df.residuals = nrx - ncx

    if (ncy == 1){
      coefficients = drop(coefficients)
      names(coefficients) = name
      fitted.values = drop(fitted.values)
      residuals = as.vector(residuals)
      effects = drop(effects)
      names(effects) = c(name, rep("", nrx-ncx))
    }else{
      rownames(effects) = c(name, rep("", nrx-ncx))
      rownames(coefficients) = name
    }

    z$add.intercept = add.intercept
    z$coefficients = coefficients
    z$residuals = residuals
    z$effects = effects
    z$rank = qr(x)$rank
    z$fitted.values = fitted.values
    z$assign = NULL
    z$qr = qr
    z$df.residual = df.residuals
    z$method = method
  }
  z
}