Modify this R code to do the following:

Create a new section after we predict final stores
In this section, we're going to create a new model function. This model function will not delete variables with high VIF. Instead, it'll automatically figure out what variables it should try to combine based on their correlation, including the original variables, the custom ones, and the per capita ones.
It'll need to figure out how to get all the variable names from both datasets even if they are of different sizes. It needs to run totally automatically when the code is run. If making it an iterative process is better (printing the summary of each new model), we can do that. Keep it all in it's own section 

Then create a new section under that new model section. This section will go through the original store prediction process but for the new store model

You can just give me the new sections that go at the bottom of the code instead of returning the full code.

Code:
# Load required packages
library(car)      # For VIF
library(dplyr)    # For data manipulation
library(corrplot) # For correlation plot
library(moments)  # For skewness calculation

# Load the data
load("C:/Users/kylek/OneDrive/Desktop/Stats 2/CroqPain_new.rda")
load("C:/Users/kylek/OneDrive/Desktop/Stats 2/CroqPainFix.rda")

# Data Preparation
CroqPainFix <- CroqPaiFix
CroqPain_new <- CroqPain_new

# Convert STOR in CroqPainFix to character (if needed)
CroqPainFix$STOR <- as.character(CroqPainFix$STOR)

# Convert EMPL in CroqPain_new to integer
CroqPain_new$EMPL <- as.integer(as.character(CroqPain_new$EMPL))

# Handle NA values in EMPL in CroqPain_new by replacing with mean
CroqPain_new$EMPL[is.na(CroqPain_new$EMPL)] <- mean(CroqPain_new$EMPL, na.rm = TRUE)

# Function to create all variables (both original and per capita)
create_variables <- function(df) {
  # Create per capita variables (excluding EARN per capita)
  df$P15_per_capita <- df$P15 / df$total
  df$P25_per_capita <- df$P25 / df$total
  df$P35_per_capita <- df$P35 / df$total
  df$P45_per_capita <- df$P45 / df$total
  df$P55_per_capita <- df$P55 / df$total
  df$COMP_per_capita <- df$COMP / df$total
  df$NCOMP_per_capita <- df$NCOMP / df$total
  df$NREST_per_capita <- df$NREST / df$total
  
  # 1. Aggregate Sales Variables
  df$Total_P_Sales <- df$P15 + df$P25 + df$P35 + df$P45 + df$P55
  
  # Percentage contribution of each price point
  df$P15_Percentage <- df$P15 / df$Total_P_Sales
  df$P25_Percentage <- df$P25 / df$Total_P_Sales
  df$P35_Percentage <- df$P35 / df$Total_P_Sales
  df$P45_Percentage <- df$P45 / df$Total_P_Sales
  
  # 2. Interaction Terms
  df$SIZE_PRICE <- df$SIZE * df$PRICE
  
  # 3. Employee Density
  df$EMPL_Density <- df$EMPL / df$SIZE
  
  # 4. Competition Intensity
  df$Competition_Density <- (df$COMP + df$NCOMP) / df$NREST
  
  # 5. Price Sensitivity
  df$Price_to_Income <- df$PRICE / df$INC
  
  # 6. Operational Efficiency Metrics
  df$Sales_per_EMPL <- df$total / df$EMPL
  
  # 7. Customer Behavior Indicators
  df$Avg_Sale_per_CLI <- df$total / df$CLI
  
  # 8. Market Penetration Metrics
  df$CLI_per_SIZE <- df$CLI / df$SIZE
  
  # 9. Log Transformations
  df$log_total <- log(df$total)
  
  # 10. Polynomial Terms
  df$SIZE_sq <- df$SIZE^2
  df$INC_sq <- df$INC^2
  
  return(df)
}

# Apply the function to both datasets
CroqPainFix <- create_variables(CroqPainFix)
CroqPain_new <- create_variables(CroqPain_new)

# Create correlation matrix for per capita variables
per_capita_vars <- CroqPainFix[, c(
  "P15_per_capita", "P25_per_capita", "P35_per_capita", 
  "P45_per_capita", "P55_per_capita", "COMP_per_capita", 
  "NCOMP_per_capita", "NREST_per_capita"
)]

# Calculate correlation matrix
cor_matrix <- cor(per_capita_vars, use = "complete.obs")
cor_matrix

# Create correlation plot
png("correlation_plot.png", width = 800, height = 800)
corrplot(cor_matrix, method = "color", type = "upper", 
         tl.col = "black", tl.srt = 45,
         addCoef.col = "black")
dev.off()

# Calculate and print skewness for per capita variables
print("Skewness Analysis for Per Capita Variables:")
skewness_results <- sapply(per_capita_vars, skewness)
for(var in names(skewness_results)) {
  skew_value <- skewness_results[var]
  is_skewed <- abs(skew_value) > 1
  skew_direction <- if(skew_value > 0) "positively" else "negatively"
  
  cat(sprintf("%s:\n", var))
  cat(sprintf("  Skewness value: %.3f\n", skew_value))
  cat(sprintf("  Status: %s\n", 
              if(is_skewed) paste("Significantly", skew_direction, "skewed") 
              else "Not significantly skewed"))
  cat("\n")
}



# Initial formula with all variables
initial_formula_custom <- EARN ~ SIZE + EMPL + INC + COMP + NCOMP + NREST + PRICE + CLI +
  Total_P_Sales + P15_Percentage + P25_Percentage + P35_Percentage + P45_Percentage +
  SIZE_PRICE + EMPL_Density + Competition_Density + Price_to_Income +
  Sales_per_EMPL + Avg_Sale_per_CLI + CLI_per_SIZE + log_total + SIZE_sq + INC_sq +
  P15_per_capita + P25_per_capita + P35_per_capita + P45_per_capita + P55_per_capita +
  COMP_per_capita + NCOMP_per_capita + NREST_per_capita

# Fit the initial model with custom variables
model_initial_custom <- lm(initial_formula_custom, data = CroqPainFix)

# Check summary of the initial model
summary(model_initial_custom)

# Check for multicollinearity using VIF
vif_values_custom <- vif(model_initial_custom)
print("VIF Values:")
print(vif_values_custom)

# Identify variables with VIF > 5
high_vif_custom <- vif_values_custom[vif_values_custom > 5]
print("Variables with VIF > 5:")
print(high_vif_custom)

# Function to iteratively remove variables with high VIF
remove_high_vif <- function(model, threshold = 5) {
  vif_values <- vif(model)
  max_vif <- max(vif_values)
  while (max_vif > threshold) {
    var_to_remove <- names(which.max(vif_values))
    cat("Removing variable due to high VIF:", var_to_remove, "\n")
    formula_terms <- attr(terms(model), "term.labels")
    formula_terms <- setdiff(formula_terms, var_to_remove)
    new_formula <- as.formula(paste("EARN ~", paste(formula_terms, collapse = " + ")))
    model <- lm(new_formula, data = CroqPainFix)
    vif_values <- vif(model)
    max_vif <- max(vif_values)
  }
  return(model)
}

# Apply the function to remove high VIF variables
model_vif_filtered_custom <- remove_high_vif(model_initial_custom, threshold = 5)

# Check summary of the model after VIF filtering
summary(model_vif_filtered_custom)

# Check VIFs after filtering
vif_values_filtered_custom <- vif(model_vif_filtered_custom)
print("VIF Values after filtering:")
print(vif_values_filtered_custom)

# Remove insignificant variables based on p-values
coeff_summary_custom <- summary(model_vif_filtered_custom)$coefficients
high_pval_vars_custom <- rownames(coeff_summary_custom)[coeff_summary_custom[, "Pr(>|t|)"] > 0.1]
high_pval_vars_custom <- setdiff(high_pval_vars_custom, "(Intercept)")
print("Variables with p-value > 0.1 to be considered for removal:")
print(high_pval_vars_custom)

# Remove variables with high p-values
final_terms_custom <- attr(terms(model_vif_filtered_custom), "term.labels")
final_terms_custom <- setdiff(final_terms_custom, high_pval_vars_custom)
final_formula_custom <- as.formula(paste("EARN ~", paste(final_terms_custom, collapse = " + ")))

# Fit the final model
model_final_custom <- lm(final_formula_custom, data = CroqPainFix)

# Summary of the final model
summary(model_final_custom)

# Check VIFs of the final model
vif_final_custom <- vif(model_final_custom)
print("Final model VIF values:")
print(vif_final_custom)

# Debugging steps to identify the source of NaN values

# More robust data preparation and cleaning
prepare_data <- function(df, is_training = TRUE) {
  # 1. Handle EMPL specifically first
  if("EMPL" %in% names(df)) {
    # Convert EMPL to numeric, keeping original values
    df$EMPL <- as.numeric(as.character(df$EMPL))
    
    if(is_training) {
      # Store mean EMPL from training data for later use
      assign("mean_empl", mean(df$EMPL, na.rm = TRUE), envir = .GlobalEnv)
    } else {
      # Use stored mean from training data
      if(exists("mean_empl")) {
        df$EMPL[is.na(df$EMPL)] <- mean_empl
      }
    }
  }
  
  # 2. Create derived variables
  df <- create_variables(df)  # Your existing create_variables function
  
  # 3. Handle infinite and NA values for all numeric columns
  numeric_cols <- sapply(df, is.numeric)
  for(col in names(df)[numeric_cols]) {
    # Replace infinite values with NA
    df[[col]][is.infinite(df[[col]])] <- NA
    
    # Handle NA values
    if(any(is.na(df[[col]]))) {
      if(is_training) {
        # Store mean for each column from training data
        assign(paste0("mean_", col), mean(df[[col]], na.rm = TRUE), envir = .GlobalEnv)
        df[[col]][is.na(df[[col]])] <- get(paste0("mean_", col))
      } else {
        # Use stored means from training data
        if(exists(paste0("mean_", col))) {
          df[[col]][is.na(df[[col]])] <- get(paste0("mean_", col))
        }
      }
    }
  }
  
  return(df)
}

# Apply the new preparation function to both datasets
CroqPainFix <- prepare_data(CroqPainFix, is_training = TRUE)
CroqPain_new <- prepare_data(CroqPain_new, is_training = FALSE)

# Refit the model with cleaned data
model_final_custom <- lm(final_formula_custom, data = CroqPainFix)

# Make predictions
predictions_new_custom <- predict(model_final_custom, newdata = CroqPain_new)

# Calculate performance ratios
CroqPain_new$Predicted_EARN <- predictions_new_custom
CroqPain_new$Performance_Ratio <- CroqPain_new$Predicted_EARN / CroqPain_new$K * 100

# Print results
print("Predictions after robust cleaning:")
print(CroqPain_new[, c("STOR", "Predicted_EARN", "K", "Performance_Ratio")])

# Verify no NaN values in predictions
if(any(is.na(predictions_new_custom))) {
  cat("\nWarning: NaN values still present in predictions\n")
  # Print problematic rows
  problem_rows <- which(is.na(predictions_new_custom))
  print("Problematic rows:")
  print(CroqPain_new[problem_rows, ])
}

# Validation

# Split the data into CroqPainFix_50 and CroqPainFix_10
CroqPainFix$STOR <- as.integer(CroqPainFix$STOR)
CroqPainFix_50 <- subset(CroqPainFix, STOR <= 50)
CroqPainFix_10 <- subset(CroqPainFix, STOR > 50)

# Re-estimate the final model using CroqPainFix_50
model_final_50_custom <- lm(formula(model_final_custom), data = CroqPainFix_50)

# Predict the performances of the 10 stores in CroqPainFix_10
predictions_10_custom <- predict(model_final_50_custom, newdata = CroqPainFix_10)

# Add predictions to CroqPainFix_10
CroqPainFix_10$Predicted_EARN <- predictions_10_custom
CroqPainFix_10$Performance_Ratio <- CroqPainFix_10$Predicted_EARN / CroqPainFix_10$K * 100

# View predicted performance ratios for validation stores
print("Predicted Performance Ratios for the 10 stores:")
print(CroqPainFix_10[, c("STOR", "Predicted_EARN", "K", "Performance_Ratio")])

# Determine which stores meet the performance ratio target of 26%
stores_to_open_custom <- subset(CroqPainFix_10, Performance_Ratio >= 26)
print("Stores to consider opening based on predicted performance ratio >= 26%:")
print(stores_to_open_custom[, c("STOR", "Performance_Ratio")])

# Final Prediction

# Convert STOR to character if needed in CroqPain_new
CroqPain_new$STOR <- as.character(CroqPain_new$STOR)

# Handle NA values in EMPL in CroqPain_new
CroqPain_new$EMPL <- as.integer(as.character(CroqPain_new$EMPL))
CroqPain_new$EMPL[is.na(CroqPain_new$EMPL)] <- mean(CroqPain_new$EMPL, na.rm = TRUE)

# Ensure all required variables are present in CroqPain_new
CroqPain_new <- create_variables(CroqPain_new)

# Get the names of variables used in the final model
model_vars <- all.vars(formula(model_final_custom))[-1]  # Remove response variable (EARN)

# Check if all required variables are present
missing_vars <- setdiff(model_vars, names(CroqPain_new))
if(length(missing_vars) > 0) {
  print("Missing variables in prediction dataset:")
  print(missing_vars)
  stop("Some variables are missing in the prediction dataset")
}

# Make predictions
tryCatch({
  predictions_new_custom <- predict(model_final_custom, newdata = CroqPain_new)
  
  # Add predictions to CroqPain_new
  CroqPain_new$Predicted_EARN <- predictions_new_custom
  CroqPain_new$Performance_Ratio <- CroqPain_new$Predicted_EARN / CroqPain_new$K * 100
  
  # View predicted performance ratios for new locations
  print("Predicted Performance Ratios for the new locations:")
  print(CroqPain_new[, c("STOR", "Predicted_EARN", "K", "Performance_Ratio")])
  
  # Recommend which locations to open based on performance ratio target of 26%
  locations_to_open_custom <- subset(CroqPain_new, Performance_Ratio >= 26)
  print("Recommended locations to open (Performance Ratio >= 26%):")
  print(locations_to_open_custom[, c("STOR", "Performance_Ratio")])
}, error = function(e) {
  print("Error in prediction:")
  print(e)
  
  # Print the structure of both datasets for debugging
  print("Structure of model variables:")
  print(str(CroqPainFix[, model_vars]))
  print("Structure of prediction variables:")
  print(str(CroqPain_new[, model_vars]))
})