R/olsatester.R

In jepperosenborg/olsat: Testing for OLS Assumptions

#' OLS Assumptions Tester
#'
#' @param model this is a description
#'
#' @return
#' @export
#'
#' @examples
#' olsatester(lm(Sepal.Length ~ Sepal.Width + Petal.Width, data = iris))
#'
#' @import car lmtest stats mctest formattable
olsatester <- function(model) {

  #library(formattable)
  #library(car) # NCV test
  #library(lmtest) # Breusch-Pagan, Likelihood Ratio
  #library(mctest) # Farrar-Glauber Test

  fstat<-summary(model)$fstatistic
  sig_pval <- pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)[["value"]]
  sig_pval_result <- ifelse(sig_pval > 0.05, "Unacceptable", "Acceptable")

  ass_lrtest <- lmtest::lrtest(model)$`Pr(>Chisq)`[2]
  ass_lrtest_result <- ifelse(ass_lrtest > 0.05, "Unacceptable", "Acceptable")

  if (length(model$coefficients) > 2) {
    model_matrix <- model.matrix(model)[,-1]
    ass_fgtest <- sum(mctest::imcdiag(model_matrix, model$model[1])$idiags[,"Klein"])
    ass_fgtest_result <- ifelse(ass_fgtest > 0, "Unacceptable", "Acceptable")
  } else {
    ass_fgtest <- 0
    ass_fgtest_result <- "Acceptable"
  }

  ass_bptest <- lmtest::bptest(model)$p.value[["BP"]]
  ass_bptest_result <- ifelse(ass_bptest < 0.05, "Unacceptable", "Acceptable")

  ass_ncvtest <- car::ncvTest(model)$p
  ass_ncvtest_result <- ifelse(ass_ncvtest < 0.05, "Unacceptable", "Acceptable")

  if(length(model$model[1][,1]) < 5000) {
    ass_swtest <- shapiro.test(model$model[1][,1])$p.value
  } else {
    ass_swtest <- shapiro.test(sample(model$model[1][,1], 5000))$p.value
  }
  ass_swtest_result <- ifelse(ass_swtest < 0.05, "Unacceptable", "Acceptable")

  ass_kstest <- ks.test(unique(model$model[1][,1]), "pnorm", mean(model$model[1][,1]), sd(model$model[1][,1]))$p.value
  ass_kstest_result <- ifelse(ass_kstest < 0.05, "Unacceptable", "Acceptable")

  ass_dwtest <- lmtest::dwtest(model)$p.value
  ass_dwtest_result <- ifelse(ass_dwtest < 0.05, "Unacceptable", "Acceptable")

  ass_bgtest <- lmtest::bgtest(model)$p.value
  ass_bgtest_result <- ifelse(ass_bgtest < 0.05, "Unacceptable", "Acceptable")

  assumptions <- c("MLR.1 Linearity", "MLR.1 Linearity", "MLR.3 No Perfect Collinearity (Multicollinearity)", "MLR.5 Homoscedasticity", "MLR.5 Homoscedasticity", "MLR.6 Normality of Errors (Normally distributed)", "MLR.6 Normality of Errors (Normally distributed)", "Independant error terms", "Independant error terms")
  names <- c("Model F-statistics", "Likelihood Ratio", "Farrar-Glauber", "Breusch-Pagan", "Non-Constant Error Variance", "Shapiro-Wilk", "Kolmogorov-Smirnov", "Durbin-Watson", "Breusch-Godfrey")
  results <- c(sig_pval_result, ass_lrtest_result, ass_fgtest_result, ass_bptest_result, ass_ncvtest_result, ass_swtest_result, ass_kstest_result, ass_dwtest_result, ass_bgtest_result)
  values <- c(round(sig_pval, digits = 4), round(ass_lrtest, digits = 4), round(ass_fgtest, digits = 0), round(ass_bptest, digits = 4), round(ass_ncvtest, digits = 4), round(ass_swtest, digits = 4), round(ass_kstest, digits = 4), round(ass_dwtest, digits = 4), round(ass_bgtest, digits = 4))
  valtype <- c("P-value", "P-value", "# of variables", "P-value", "P-value", "P-value", "P-value", "P-value", "P-value")
  interpretation <- c("H0: Non-linearity", "H0: Non-linearity", "Variables with collinearity", "H0: Homoscedasticity", "H0: Homoscedasticity", "H0: Normally distributed", "H0: Normally distributed", "H0: Independance", "H0: Independance")
  output_frame <- data.frame("Assumption" = assumptions, "Test" = names, "Values" = values, "Value type" = valtype, "Interpretation" = interpretation, "Acceptance" = results)

  acceptance_formatter <- formatter("span",
                                    style = x ~ style(
                                      font.weight = "bold",
                                      color = ifelse(x == "Acceptable", "#71CA97", ifelse(x == "Unacceptable", "#ff7f7f", "black"))))

  output <- formattable::formattable(output_frame, align = "l", list(`Acceptance` = acceptance_formatter))


  return(output)
}