R/ROB_Statistics.R

In ZhenWei10/Sherry-Chapter1:

#'@title A function that run NLR and calculate R squared over the model matrix
#'
#'@description A function that calculate the R squared of a NLR for RV prediction.
#'
#'@param MM a matrix of the model matrix.
#'
#'@return a vector for NLR R squares with 14 different designs.
#'
#'@import MASS
#'@export
#'
ROB_Statistics <- function(MM){

  return_matrix = matrix(NA,nrow = 3, ncol = 4)
  rownames(return_matrix) = c("coef_est","pvalue","adj_Rsquared")
  colnames(return_matrix) = c("V_nof","TN_nof","V_filter","TN_filter")

  MM = MM[,colnames(MM) != "X"]
  indx_basic = grepl("Monday|Y",colnames(MM))
  MM = na.omit(MM)
  MM = as.data.frame(scale(MM))

  #define a helper function
  rlm_stat <- function(Model_M){
    res <- rlm(Y ~ ., data = Model_M)
    statistics <- summary(res)
    R2 <-  summary( lm(Model_M$Y~fitted(res) - 1) )$r.squared
    n = nrow(Model_M)
    p = ncol(Model_M) - 1
    adj_R2 = 1-(1-R2)*((n-1)/(n-p-1))
    pvalue <- 2*(1-pt(abs(statistics$coefficients[3,3]),df = n-p))
    coef_est <- statistics$coefficients[3,1]
    return(c(coef_est,pvalue,adj_R2))
  }

  return_matrix[,1] =  rlm_stat(MM[,indx_basic|colnames(MM) == "V_raw"])
  return_matrix[,2] =  rlm_stat(MM[,indx_basic|colnames(MM) == "TN_raw"])
  return_matrix[,3] =  rlm_stat(MM[,indx_basic|colnames(MM) == "V_weekly_stl"])
  return_matrix[,4] =  rlm_stat(MM[,indx_basic|colnames(MM) == "TN_weekly_stl"])

  return(return_matrix)
}