assignment4: Mutlpile Regression

#' Illustration of Linear Regression
#' @field  formula is the formula object of R , which independent and dependent variables names
#' @field  data it is a data.frame
#' @description This RC object calculates multiple regression
#' @references https://en.wikipedia.org/wiki/Linear_regression
#' @export linreg
#' @exportClass linreg
#' 
#'  



linreg <- setRefClass("linreg" , 
                      fields = list(formula = "formula", 
                                    data = "data.frame",
                                    residuals = "vector",
                                    coefficients = "vector",
                                    output_y = "vector",
                                    coef_with_names = "vector",
                                    data_frame_name = "character",
                                    summary_matrix = "matrix",
                                    df = "numeric",
                                    RSD = "numeric",
                                    parameters = "matrix"),
                      methods = list(initialize = function(formula, data){
                        
                        stopifnot(!is.null(data))
                        #stopifnot(class(formula) == "formula",NULL)          
                       
                        .self$formula <<- formula
                        .self$data <<- data
                       # browser()
                        data_frame_name <<- deparse(substitute(data))
                        .self$regression()
                        
                      },
                      regression = function(){

                          X <- model.matrix(formula, data)
                          
                          y_col <- data[all.vars(formula)[1]]
                          y <- as.matrix(y_col)
                          #browser()
                          transpose_X <- t(X)

                          X_tranp_inversion <- solve(transpose_X %*% X)
                          
                          beta_hat <- X_tranp_inversion %*% (transpose_X %*% y) # also called estimates
                          .self$coefficients <<- as.vector(beta_hat)
                          .self$parameters <<- beta_hat
                          #browser()
                          fitted_value_y_hat <- X %*% beta_hat
                          .self$output_y <<- as.vector(fitted_value_y_hat)
                          
                          resid <-  y- fitted_value_y_hat
                          .self$residuals <<- as.vector(resid)
                          
                          degree_freedom <- nrow(X) - ncol(X)
                          df <<- degree_freedom
                          variance_sqr <- (t(resid) %*% resid) / degree_freedom
                          variance <- sqrt(variance_sqr)
                          RSD <<- as.numeric(variance)
                          variance_beta <- as.vector(variance_sqr) * diag(X_tranp_inversion)
                          std_error <- sqrt(variance_beta)
                          #parameters <- diag(X_tranp_inversion)
                          #parameters_names <- names(parameters)
                          #browser()
                          t_value <- beta_hat/ sqrt(variance_beta)
                          
                          #p_value <- 2*pt(-abs(t_value),degree_freedom)
                          
                          p_value <- 2 * (1 - pt(abs(t_value),degree_freedom))
                          
                          
                          .self$coef_with_names <<- as.vector(beta_hat) 
                          names(.self$coef_with_names) <- row.names(beta_hat)
                          
                          colum_names <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
                          .self$summary_matrix <<- matrix(.self$coef_with_names)
                          .self$summary_matrix <<- cbind(.self$summary_matrix,std_error)
                          .self$summary_matrix <<- cbind(.self$summary_matrix,t_value)
                          .self$summary_matrix <<- cbind(.self$summary_matrix,p_value)
                          colnames(.self$summary_matrix) <- colum_names
                          #browser()
                          #dim(.self$summary_matrix)
                      },
                      show = function(){
                        .self$print()
                      },
                      print = function(){
                        
                        formula_vars <- as.character(formula)
                        #dataframe_name <- deparse(substitute(.self$data)) 
                        #browser()
                        cat("Call:\n")
                        cat(paste0("linreg(formula = ",
                                   formula_vars[2],
                                   " ",
                                   formula_vars[1],
                                   " ",
                                   formula_vars[3],
                                   ", data = ",
                                   .self$data_frame_name ,")"))
                        
                        col_names <- names(.self$coef_with_names)
                        
                        col_values <- .self$coef_with_names
                        names(col_values) <- NULL
                        col_values <- round(col_values,3)
                        
                        #browser()
                        cat("\n\nCoefficients:")
                        cat("\n")
                        cat("    ",col_names,sep = "  ")
                        cat("\n")
                        values_formatting = "      "
                        for (i in 1:length(col_values)) {
                          
                          values_formatting<- paste0(values_formatting,
                              format(col_values[i] ,
                                     justify = "right",
                                     width = max(nchar(col_names[i]),
                                                 nchar(col_values[i]))
                                     ))
                          values_formatting <- paste0(values_formatting,"  ")    
                      
                        }
                        cat(values_formatting)
                        #cat(paste0("    " ,col_values , sep = "  " , fill = TRUE))
                        
                      },
                      summary = function(){
                        
                        name <- "Coefficients:"
                        output_dataframe <- as.data.frame(.self$summary_matrix)
                        last_line <- (paste0("Residual standard error: ", RSD," on ",  df, " degrees of freedom"))
                        
                        
                        #summary_list <- list(name , output_dataframe , last_line)
                        #class(summary_list) <- "summary.linreg"
                        cat(paste0(name, "\n"))
                        cat("\t" , colnames(output_dataframe) , "\n" , sep = "\t")
                        for (item in 1:(nrow(output_dataframe))) {
                          #browser()
                          values <- as.vector(unlist(output_dataframe[item,]))
                          row_name <- rownames(output_dataframe)[item]
                          #browser()
                          cat(row_name , values ,collapse = " ***" , "\n")
                          #browser()
                        }
                        cat("\n" ,last_line)
                        
                      },
                      resid =  function(){
                        return(.self$residuals)
                      },
                      pred =  function(data = iris){
                        stopifnot(inherits(data,"data.frame"))
                        #stopifnot(inherits(linreg_obj,"linreg"))
                        
                        new_X <- model.matrix(.self$formula, data)
                        prediction <- new_X %*% parameters
                        return(as.vector(prediction))
                      },
                      coef = function(){
                        return(.self$coef_with_names)
                      },
                      plot = function(){
                        #x <- as.factor(.self$output_y) 
                        #y <- as.factor(.self$residuals)
                        require(ggplot2)
                        require(gridExtra)
                        
                        formula_string <- paste0(as.character(.self$formula)[2],
                                                 as.character(.self$formula)[1],
                                                 as.character(.self$formula)[3]) 
                        
                        plot_data <- data.frame(resid = .self$resid() , fitted = .self$pred())
                        plot_data["std_residuals"] = sqrt(.self$resid())
                        
                        # This is first plot for residual and fitted values
                        
                        # """ For debugging """
                        #browser()
                       
                        plot1 <-  ggplot(data = plot_data) + 
                          aes(x = fitted , y = resid) +
                          geom_point(shape=1, size=2)+
                          labs(title = "Residuals vs Fitted",
                                        x = "Fitted values",
                                        y = "Residuals",
                                        caption = paste0("linreg( ",formula_string," )"))+
                          theme(
                            plot.title = element_text(hjust = 0.5),
                            plot.caption =  element_text(hjust = 0.5),
                            plot.background = element_rect(fill = "transparent" , colour = NA),
                            panel.grid.minor = element_blank(),
                            panel.grid.major = element_blank(),
                            panel.border = element_rect(colour = "black" , fill = NA , size = 1)
                          ) +
                          # geom_hline(yintercept = 0,
                          #                     color="grey",
                          #                     linetype = "dashed",
                          #                     size = 0.3)+
                          geom_smooth(method = "loess" ,formula = y ~ x,se =FALSE)
                        
                        # This is second plot for residual and fitted values
                        
                       plot2 <-  ggplot(data = na.omit(plot_data)) +
                          aes(x = fitted , y = std_residuals ) +
                          geom_point(shape=1, size=2)+
                          labs(title = "Scaled Location",
                                        x = "Fitted values",
                                        y = expression(paste(sqrt("Standardized Residuals"))),
                                        caption = paste0("linreg( " ,formula_string," )"))+
                          theme(
                            plot.title = element_text(hjust = 0.5),
                            plot.caption =  element_text(hjust = 0.5),
                            plot.background = element_rect(fill = "transparent" , colour = NA),
                            panel.grid.minor = element_blank(),
                            panel.grid.major = element_blank(),
                            panel.border = element_rect(colour = "black" , fill = NA , size = 1)
                          ) +

                          geom_hline(yintercept = 0,
                                              color="grey",
                                              linetype = "dashed",
                                              size = 0.3)+
                          geom_smooth(method = "loess" ,formula = y ~ x,se =FALSE)
                       grid.arrange(plot1, 
                                    plot2,
                                    ncol = 1,
                                    nrow = 2)
                        }
                                     )
                                     
                      )