R/linreg.R

In harjew/lab4G3: Computer lab 4

#' Linear regression using least squares
#' 
#' linear regression using least square method.
#' 
#' @param formula An object of class "formula".
#' @param data An object of class "data frame".
#' 
#' 
#' @export



linreg<-function(formula,data){
  
  stopifnot(class(formula) == "formula",is.data.frame(data))
  # Dataframe with dummys
 
  data_for_lm<-stats::model.frame(formula,data=data)
  dummy = FALSE
  
  for(i in 1:length(data_for_lm)){
    if(class(data_for_lm[,i])== "factor"){
    dummy<-TRUE
    }
  }  
  
  

  if(dummy==TRUE){
  data_for_lm<-fastDummies::dummy_cols(data_for_lm,remove_first_dummy = T)
  k<-length(data_for_lm)
  i<-1
  while(i< k){
    if(class(data_for_lm[,i])== "factor"){
      data_for_lm[,i]<-NULL
      k<-k-1
    }
    i<-i+1
  }
  }
  
  
  # create the matrices
  y<-matrix(data_for_lm[,1],ncol = 1)
  x_v<-rep(1,length(data_for_lm[,1]))
  for(i in 2:length(data_for_lm)){
    x_v<-c(x_v,data_for_lm[,i])
  }
  x<-matrix(x_v,ncol = length(data_for_lm))
  # My lm list
  # call
  my_lm<-list()
  ca<-match.call()
  my_lm$call<-c(ca)
    
  # Beta
  Beta<-(solve(t(x)%*%x))%*%t(x)%*%y
  beta_names<-"(Intercept)"
  for(i in 2:length(data_for_lm)){
    beta_names[i]<-names(data_for_lm[i])
  }
  #print(beta_names)
  Beta_saved<-as.vector(Beta)
  names(Beta_saved)<-beta_names
  my_lm$Coefficients<-Beta_saved
  # The fitted values
  y_hat<-x%*%Beta
  y_saved<-as.vector(y_hat)
  names(y_saved)<-1:length(y_saved)
  my_lm$fitted_values<-y_saved
  

  # Residuals
  e<-y-y_hat
  e_saved<-as.vector(e)
  names(e_saved)<-1:length(e_saved)
  my_lm$Residuals<-e_saved
  
  
  # Defrees of freedom 
  df<-length(data_for_lm[,1]) - length(data_for_lm)
  my_lm$df<-df
  
  #R squared and adjusted R scuared
  sst<-(as.vector(y-mean(y)))^2
  ssr<-((y-y_saved)^2)
  r_squared<-(1-sum(ssr)/sum(sst))
  my_lm$r_squared<-r_squared
  adj_r_squared<-1-(1-r_squared)*((length(data_for_lm[,1])-1))/(df-1)
  my_lm$adj_r_squared<-adj_r_squared
  
  # Residuals variance
  sigma_2<-as.vector(t(e)%*%e/df)
  my_lm$Residuals_variance<-sigma_2

  # Variance of regressin coefficietns
  var_Beta<-sigma_2*solve(t(x)%*%x)

  var_Beta_estimate<-vector()
  
  for(i in 1:dim(var_Beta)[2]){
    var_Beta_estimate[i]<-var_Beta[i,i]
  }
  se_error<-(var_Beta_estimate^.5)
  
  my_lm$Std.Error<-(se_error)
  # t-values for each coeffissient
  
  t_Beta<-Beta/(se_error)
  my_lm$t_values<-as.vector(t_Beta)
  # p-value
  p_value<-vector()
  for (i in 1:length(t_Beta)){
    p_value[i]<-2*stats::pt(-abs(t_Beta[i]),df=df)
  }
  my_lm$p_value<-p_value
  
  LM_Shape <- function(data){
    base::structure(data, class = "linreg")
  }
  return(LM_Shape(my_lm))
  
}