r R中的自定义评估回归模型脚本用于Azure机器学习

Posted 2021-05-24

# Map 1-based optional input ports to variables
df <- maml.mapInputPort(1) # class: data.frame
df_orig <- maml.mapInputPort(2)

### INPUT PARAMETERS #################################################
# k is the number of parameters to be estimated excuding the intercept
target_variable_name            <- "LOG_QUANTITY"
predicted_values_variable_name  <- "Scored Labels"
k                               <- 23
######################################################################


### LIBRARIES ##################################
library(ggplot2)


### FUNCTIONS ##################################

# QQplot as ggplot
# Source: https://gist.github.com/rentrop/d39a8406ad8af2a1066c
gg_qq <- function(x, distribution = "norm", ..., line.estimate = NULL, conf = 0.95,
                  labels = names(x)){
  q.function <- eval(parse(text = paste0("q", distribution)))
  d.function <- eval(parse(text = paste0("d", distribution)))
  x <- na.omit(x)
  ord <- order(x)
  n <- length(x)
  P <- ppoints(length(x))
  df <- data.frame(ord.x = x[ord], z = q.function(P, ...))

  if(is.null(line.estimate)){
    Q.x <- quantile(df$ord.x, c(0.25, 0.75))
    Q.z <- q.function(c(0.25, 0.75), ...)
    b <- diff(Q.x)/diff(Q.z)
    coef <- c(Q.x[1] - b * Q.z[1], b)
  } else {
    coef <- coef(line.estimate(ord.x ~ z))
  }

  zz <- qnorm(1 - (1 - conf)/2)
  SE <- (coef[2]/d.function(df$z)) * sqrt(P * (1 - P)/n)
  fit.value <- coef[1] + coef[2] * df$z
  df$upper <- fit.value + zz * SE
  df$lower <- fit.value - zz * SE

  if(!is.null(labels)){ 
    df$label <- ifelse(df$ord.x > df$upper | df$ord.x < df$lower, labels[ord],"")
    }

  p <- ggplot(df, aes(x=z, y=ord.x)) +
    geom_point(color="blue") + 
    geom_abline(intercept = coef[1], slope = coef[2], linetype = "dashed") +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2) 
  if(!is.null(labels)) p <- p + geom_text( aes(label = label))
  
  #print(p)
  #coef
  return(list(plot = p, coef))
}



### MAIN #######################################

#---------------
#--- Metrics ---
#---------------
actual_values <- df[, target_variable_name]
predicted_values <- df[, predicted_values_variable_name]

error_values <- actual_values - predicted_values

sd_error_values <- sd(error_values)
standardized_residual_values <- error_values / sd_error_values

n <- length(actual_values)
degrees_of_freedom <- n - k - 1

# Names of the resulting data frame's rows
names <- c("MAE", "MAPE", "RMSE", "NRMSE", "R^2", "adj-R^2", "S")


mae     <- c(value = mean( abs(error_values) ))
mape    <- c(value = mean( abs(100 * error_values / actual_values) ))
rmses   <- c(value = sqrt( mean( error_values^2 , na.rm = TRUE ) ))
nrmses  <- c(value = rmses /
                ( max(actual_values, na.rm = TRUE) - min(actual_values, na.rm = TRUE) )
           )
rsq     <- c( value = 1 - ( sum(error_values[!is.na(error_values)]^2) /
                            sum( (actual_values[!is.na(actual_values)] - mean(actual_values, na.rm = TRUE) )^2) )
           )
s       <- c(value = sqrt( sum(error_values[!is.na(error_values)]^2) / degrees_of_freedom))
radjsq  <- 1 - (1 - rsq) * (n - 1) / (n - k - 1)

res <- as.data.frame(
        cbind(stats = names,
              round(rbind(mae, mape, rmses, nrmses, rsq, radjsq, s), 6)
        )
       )


#-------------
#--- Plots ---
#-------------

title0 <- 'HISTOGRAM OF RESIDUALS'
subtitle0 <- target_variable_name

p0 <- ggplot(data = df, aes(error_values)) +
        geom_histogram(fill = "blue") +
        ggtitle(bquote(atop(.(title0), atop(italic(.(subtitle0)), "")))) +
        xlab("Residuals") + ylab("Frequence")

#print(p0)  # Use to debug the code
try(print(p0), silent = TRUE)



title1 <- 'PREDICTED VS ACTUAL'
subtitle1 <- target_variable_name

p1 <- ggplot(data = df, aes(x = predicted_values, y = actual_values)) +
        geom_point(color="blue") +
        geom_abline(intercept = 0, linetype = "dashed") +
        ggtitle(bquote(atop(.(title1), atop(italic(.(subtitle1)), "")))) +
        xlab("Predicted values") + ylab("Actual values") #+ expand_limits(x = 0, y = 0)

#print(p1)  # Use to debug the code
try(print(p1), silent = TRUE)



title2 <- 'PREDICTED VS RESIDUAL'
subtitle2 <- target_variable_name

p2 <- ggplot(data = df, aes(x = predicted_values, y = standardized_residual_values)) +
        geom_point(color="blue") +
        geom_hline(yintercept = 0, linetype = "dashed") +
        ggtitle(bquote(atop(.(title2), atop(italic(.(subtitle2)), "")))) +
        xlab("Predicted values") + ylab("Standardized residuals")

#print(p2)  # Use to debug the code
try(print(p2), silent = TRUE)



title3 <- 'NORMAL Q-Q PLOT'
subtitle3 <- target_variable_name

p3 <- gg_qq(error_values)$plot +
        ggtitle(bquote(atop(.(title3), atop(italic(.(subtitle3)), "")))) +
        xlab("Theoretical Quantiles") + ylab("Residuals Quantiles")

#print(p3)  # Use to debug the code
try(print(p3), silent = TRUE)




for (regressor_name in names(df_orig))
{
  if (regressor_name %in% names(df))
    feature_type <- "FEATURE"
  else
    feature_type <- "EXTERNAL FEATURE"
    
  variable_nature <- ifelse(regressor_name == target_variable_name, "LABEL", feature_type)
  if (!grepl("Scored", regressor_name, fixed = TRUE))
  {
    title4 <- paste(variable_nature, ' VS RESIDUAL', sep="")
    subtitle4 <- regressor_name
    
    if (class(df_orig[,regressor_name]) != "factor")
    {
      p4 <- ggplot(data = df_orig, aes(x = get(regressor_name), y = standardized_residual_values)) +
              geom_point(color="blue") +
              geom_hline(yintercept = 0, linetype = "dashed") +
              ggtitle(bquote(atop(.(title4), atop(italic(.(subtitle4)), "")))) +
              xlab(tolower(regressor_name)) + ylab("Standardized residuals")
    }
    else
    {
      p4 <- ggplot(data = df_orig, aes(x = get(regressor_name), y = standardized_residual_values)) +
              geom_boxplot(color="blue") +
              #geom_hline(yintercept = 0, linetype = "dashed") +
              ggtitle(bquote(atop(.(title4), atop(italic(.(subtitle4)), "")))) +
              xlab(tolower(regressor_name)) + ylab("Standardized residuals")
    }
    print(p4)  # Use to debug the code
    #try(print(p4), silent = TRUE)
  }
}

# Select data.frame to be sent to the output Dataset port
maml.mapOutputPort("res");