R/linear_model.R

In damingli09/bis557: Package for homework of BIS557 Computational Statistics

#' @title Fits linear model
#' @description Fits linear model by solving normal equation
#' @param form the input formula object
#' @param dat the input dataframe
#' @param constrasts list of contrasts for factor input
#' @example 
#' fit <- linear_model(Sepal.Length ~ ., iris, contrasts = list(Species = "contr.sum"))
#' print(fit)
#' @export

linear_model <- function(form, dat, contrasts=NULL){
    if (! is.null(contrasts)){
      contrasts(dat[,names(contrasts)]) <- contrasts[[names(contrasts)]]
    }
  
    X <- model.matrix(form,dat)
    y_name <- as.character(form)[2]
    y <- matrix(model.frame(form,dat)[,y_name], ncol=1)
    
    #betas <- solve(t(X)%*%X)%*%t(X)%*%y
    betas <- solve.qr(qr(X), y)
    betas[which(betas==0)] = NA
    betas_name <- rownames(betas)
    betas <- as.numeric(betas)
    names(betas) <- betas_name
    
    res <- list(coefficients = betas)
    #class(res) <- "my_lm"
    res
}