This was my result, any ideas on what custom variables might be good to add?

Stats Project Redo
Kyle Kochanowski
2024-10-04


Libraries and data
# Load required packages
library(car)      # For VIF
## Loading required package: carData
library(dplyr)    # For data manipulation
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(corrplot) # For correlation plot
## corrplot 0.94 loaded
library(leaps)    # For variable selection (if needed)

# 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 prep
#Ensure data frames are correctly named
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)


Building models
# Step 1: Explore the data
# View structure of CroqPainFix
str(CroqPainFix)
## tibble [65 × 17] (S3: tbl_df/tbl/data.frame)
##  $ STOR : chr [1:65] "1" "2" "3" "4" ...
##  $ EARN : num [1:65] 28.34 -1.46 68.93 202.07 115.78 ...
##  $ K    : num [1:65] 861 630 1074 882 931 ...
##  $ SIZE : int [1:65] 129 91 140 184 144 160 94 100 85 92 ...
##  $ EMPL : int [1:65] 14 12 13 7 14 11 5 8 14 6 ...
##  $ total: int [1:65] 8580 8460 19250 20920 11660 25780 19000 18500 14210 17440 ...
##  $ P15  : int [1:65] 980 1290 2940 3570 1700 4640 3600 3450 1930 3520 ...
##  $ P25  : int [1:65] 1280 720 2490 4930 1140 3150 2330 2560 4280 1780 ...
##  $ P35  : int [1:65] 560 1200 3710 4420 2200 5720 4750 3630 1740 4350 ...
##  $ P45  : int [1:65] 1000 1490 4030 4300 2140 5330 4970 3520 2060 4020 ...
##  $ P55  : int [1:65] 3100 3100 5270 2960 2630 5920 3030 4800 2960 3470 ...
##  $ INC  : num [1:65] 27.6 28.3 30.2 27.6 33.9 32.5 33.1 29.7 28.4 28.3 ...
##  $ COMP : int [1:65] 8 2 5 2 1 3 0 3 4 8 ...
##  $ NCOMP: int [1:65] 1 2 4 1 3 9 11 3 12 1 ...
##  $ NREST: int [1:65] 45 27 5 7 25 8 89 14 43 76 ...
##  $ PRICE: num [1:65] 16.1 11.4 21.7 11.8 16.6 22.1 24.3 16.4 12.9 13 ...
##  $ CLI  : int [1:65] 129 116 142 138 126 137 134 132 129 117 ...
##  - attr(*, "description")= chr "## Croq'Pain\n\nData on store earnings and other characteristics for the Croq'Pain case. There are total of 60 "| __truncated__
# Summary statistics
summary(CroqPainFix)
##      STOR                EARN              K               SIZE      
##  Length:65          Min.   :-42.98   Min.   : 504.3   Min.   : 52.0  
##  Class :character   1st Qu.: 77.56   1st Qu.: 647.8   1st Qu.: 83.0  
##  Mode  :character   Median :140.50   Median : 761.3   Median : 95.0  
##                     Mean   :142.45   Mean   : 836.6   Mean   :115.8  
##                     3rd Qu.:202.07   3rd Qu.: 906.5   3rd Qu.:141.0  
##                     Max.   :399.24   Max.   :1719.5   Max.   :285.0  
##       EMPL           total            P15            P25            P35      
##  Min.   : 4.00   Min.   : 6040   Min.   : 720   Min.   :  60   Min.   : 250  
##  1st Qu.: 8.00   1st Qu.:13550   1st Qu.:2050   1st Qu.:1800   1st Qu.:2360  
##  Median :11.00   Median :19920   Median :3310   Median :2810   Median :3730  
##  Mean   :10.88   Mean   :18891   Mean   :3104   Mean   :2804   Mean   :3710  
##  3rd Qu.:14.00   3rd Qu.:23400   3rd Qu.:4050   3rd Qu.:3820   3rd Qu.:5190  
##  Max.   :17.00   Max.   :30020   Max.   :5150   Max.   :5540   Max.   :6280  
##       P45            P55            INC             COMP       
##  Min.   : 580   Min.   : 910   Min.   :25.20   Min.   : 0.000  
##  1st Qu.:2320   1st Qu.:3440   1st Qu.:28.70   1st Qu.: 1.000  
##  Median :3970   Median :4450   Median :30.00   Median : 2.000  
##  Mean   :3727   Mean   :4418   Mean   :31.18   Mean   : 3.138  
##  3rd Qu.:5110   3rd Qu.:5670   3rd Qu.:33.80   3rd Qu.: 5.000  
##  Max.   :5990   Max.   :6720   Max.   :38.80   Max.   :10.000  
##      NCOMP            NREST           PRICE            CLI       
##  Min.   : 0.000   Min.   : 5.00   Min.   : 7.60   Min.   :110.0  
##  1st Qu.: 3.000   1st Qu.:13.00   1st Qu.:12.30   1st Qu.:122.0  
##  Median : 5.000   Median :25.00   Median :15.90   Median :129.0  
##  Mean   : 6.092   Mean   :29.83   Mean   :16.89   Mean   :128.1  
##  3rd Qu.: 9.000   3rd Qu.:44.00   3rd Qu.:20.60   3rd Qu.:134.0  
##  Max.   :14.000   Max.   :96.00   Max.   :33.40   Max.   :142.0
# Check for missing values
colSums(is.na(CroqPainFix))
##  STOR  EARN     K  SIZE  EMPL total   P15   P25   P35   P45   P55   INC  COMP 
##     0     0     0     0     0     0     0     0     0     0     0     0     0 
## NCOMP NREST PRICE   CLI 
##     0     0     0     0
# Step 2: Examine correlations among variables
# Select numeric variables for correlation analysis
numeric_vars <- sapply(CroqPainFix, is.numeric)
cor_matrix <- cor(CroqPainFix[, numeric_vars], use = "complete.obs")

# Visualize correlation matrix
corrplot(cor_matrix, method = "circle", tl.cex = 0.7)


# Step 3: Build initial regression model with all potential explanatory variables
# Define the formula with all variables (adjust variable names as per your dataset)
initial_formula <- EARN ~ SIZE + EMPL + total + P15 + P25 + P35 + P45 + P55 +
  INC + COMP + NCOMP + NREST + PRICE + CLI

# Fit the initial model
model_initial <- lm(initial_formula, data = CroqPainFix)

# Check summary of the initial model
summary(model_initial)
## 
## Call:
## lm(formula = initial_formula, data = CroqPainFix)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -71.042 -20.590   5.543  22.047  68.578 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.724e+02  9.041e+01  -4.119 0.000143 ***
## SIZE         7.736e-01  9.443e-02   8.192 8.41e-11 ***
## EMPL        -1.138e+00  1.397e+00  -0.814 0.419425    
## total       -1.083e-02  1.196e-02  -0.906 0.369378    
## P15          5.756e-02  2.647e-02   2.175 0.034403 *  
## P25          1.393e-02  1.260e-02   1.105 0.274296    
## P35          1.423e-02  2.105e-02   0.676 0.502038    
## P45          2.306e-03  3.019e-02   0.076 0.939402    
## P55          1.077e-02  1.333e-02   0.808 0.422787    
## INC          8.811e+00  1.584e+00   5.562 1.04e-06 ***
## COMP        -3.090e+00  2.218e+00  -1.393 0.169787    
## NCOMP       -7.669e-01  1.468e+00  -0.523 0.603580    
## NREST        1.475e+00  2.313e-01   6.379 5.61e-08 ***
## PRICE       -3.110e+00  8.931e-01  -3.483 0.001041 ** 
## CLI          4.974e-01  6.348e-01   0.784 0.436936    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.93 on 50 degrees of freedom
## Multiple R-squared:  0.8631, Adjusted R-squared:  0.8248 
## F-statistic: 22.52 on 14 and 50 DF,  p-value: < 2.2e-16
# Step 4: Check for multicollinearity using VIF
vif_values <- vif(model_initial)
print("VIF Values:")
## [1] "VIF Values:"
print(vif_values)
##       SIZE       EMPL      total        P15        P25        P35        P45 
##   1.152858   1.124262 251.351159  45.670193  12.292057  60.500390 103.262731 
##        P55        INC       COMP      NCOMP      NREST      PRICE        CLI 
##  16.327545   1.183307   1.305389   1.215359   1.217932   1.388799   1.323541
# Identify variables with VIF > 5
high_vif <- vif_values[vif_values > 5]
print("Variables with VIF > 5:")
## [1] "Variables with VIF > 5:"
print(high_vif)
##     total       P15       P25       P35       P45       P55 
## 251.35116  45.67019  12.29206  60.50039 103.26273  16.32754
# Step 5: Iteratively remove variables with high VIF
# Function to iteratively remove variables with highest 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")
    # Update the formula
    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 = " + ")))
    # Refit the model
    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 <- remove_high_vif(model_initial, threshold = 5)
## Removing variable due to high VIF: total 
## Removing variable due to high VIF: P35 
## Removing variable due to high VIF: P15
# Check summary of the model after VIF filtering
summary(model_vif_filtered)
## 
## Call:
## lm(formula = new_formula, data = CroqPainFix)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -68.871 -25.795   3.155  22.789  77.503 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.817e+02  9.242e+01  -4.130 0.000129 ***
## SIZE         7.220e-01  9.509e-02   7.593 5.01e-10 ***
## EMPL        -1.095e+00  1.436e+00  -0.762 0.449170    
## P25          4.439e-03  4.228e-03   1.050 0.298626    
## P45          2.723e-02  4.701e-03   5.793 3.86e-07 ***
## P55          4.212e-03  4.840e-03   0.870 0.388066    
## INC          9.689e+00  1.591e+00   6.089 1.31e-07 ***
## COMP        -2.652e+00  2.281e+00  -1.163 0.250118    
## NCOMP       -1.119e+00  1.475e+00  -0.759 0.451245    
## NREST        1.440e+00  2.357e-01   6.111 1.21e-07 ***
## PRICE       -3.430e+00  9.071e-01  -3.781 0.000398 ***
## CLI          3.734e-01  6.502e-01   0.574 0.568232    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.12 on 53 degrees of freedom
## Multiple R-squared:  0.8454, Adjusted R-squared:  0.8133 
## F-statistic: 26.34 on 11 and 53 DF,  p-value: < 2.2e-16
# Check VIFs after filtering
vif_values_filtered <- vif(model_vif_filtered)
print("VIF Values after filtering:")
## [1] "VIF Values after filtering:"
print(vif_values_filtered)
##     SIZE     EMPL      P25      P45      P55      INC     COMP    NCOMP 
## 1.096746 1.113671 1.298269 2.349379 2.020361 1.119824 1.295011 1.151215 
##    NREST    PRICE      CLI 
## 1.186942 1.344258 1.302850
# Step 6: Remove insignificant variables based on p-values
# Identify variables with p-value > 0.1
coeff_summary <- summary(model_vif_filtered)$coefficients
high_pval_vars <- rownames(coeff_summary)[coeff_summary[, "Pr(>|t|)"] > 0.1]
# Exclude Intercept
high_pval_vars <- setdiff(high_pval_vars, "(Intercept)")
print("Variables with p-value > 0.1 to be considered for removal:")
## [1] "Variables with p-value > 0.1 to be considered for removal:"
print(high_pval_vars)
## [1] "EMPL"  "P25"   "P55"   "COMP"  "NCOMP" "CLI"
# Remove variables with high p-values
final_terms <- attr(terms(model_vif_filtered), "term.labels")
final_terms <- setdiff(final_terms, high_pval_vars)
final_formula <- as.formula(paste("EARN ~", paste(final_terms, collapse = " + ")))

# Fit the final model
model_final <- lm(final_formula, data = CroqPainFix)

# Summary of the final model
summary(model_final)
## 
## Call:
## lm(formula = final_formula, data = CroqPainFix)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.739 -32.101   4.533  27.260  69.348 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.495e+02  5.030e+01  -6.948 3.31e-09 ***
## SIZE         7.461e-01  9.260e-02   8.057 4.38e-11 ***
## P45          3.227e-02  3.099e-03  10.413 5.52e-15 ***
## INC          9.256e+00  1.537e+00   6.023 1.19e-07 ***
## NREST        1.382e+00  2.176e-01   6.351 3.36e-08 ***
## PRICE       -2.641e+00  7.800e-01  -3.386  0.00127 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37.87 on 59 degrees of freedom
## Multiple R-squared:  0.8302, Adjusted R-squared:  0.8158 
## F-statistic: 57.68 on 5 and 59 DF,  p-value: < 2.2e-16
# Check VIFs of the final model
vif_final <- vif(model_final)
print("Final model VIF values:")
## [1] "Final model VIF values:"
print(vif_final)
##     SIZE      P45      INC    NREST    PRICE 
## 1.054300 1.034647 1.058910 1.025606 1.007437


Validation
# Split the data into CroqPainFix_50 and CroqPainFix_10
CroqPainFix$STOR <- as.integer(CroqPainFix$STOR)  # Ensure STOR is integer for comparison

CroqPainFix_50 <- subset(CroqPainFix, STOR <= 50)
CroqPainFix_10 <- subset(CroqPainFix, STOR > 50)

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

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

# Add predictions to CroqPainFix_10
CroqPainFix_10$Predicted_EARN <- predictions_10

# Calculate predicted performance ratio
# Assuming performance ratio is EARN / K (investment)
CroqPainFix_10$Performance_Ratio <- CroqPainFix_10$Predicted_EARN / CroqPainFix_10$K * 100  # in percentage

# View predicted performance ratios
print("Predicted Performance Ratios for the 10 stores:")
## [1] "Predicted Performance Ratios for the 10 stores:"
print(CroqPainFix_10[, c("STOR", "Predicted_EARN", "K", "Performance_Ratio")])
## # A tibble: 10 × 4
##     STOR Predicted_EARN     K Performance_Ratio
##    <int>          <dbl> <dbl>             <dbl>
##  1    51          151.   776.             19.5 
##  2    52          157.   648.             24.2 
##  3    53           60.3  690.              8.75
##  4    54          105.   715.             14.7 
##  5    55           54.4  650.              8.37
##  6    56           53.0  788.              6.72
##  7    57          236.   782.             30.2 
##  8    58          260.  1558.             16.7 
##  9    59          131.   936.             14.0 
## 10    60          253.   688.             36.7
# Determine which stores meet the performance ratio target of 26%
stores_to_open <- subset(CroqPainFix_10, Performance_Ratio >= 26)
print("Stores to consider opening based on predicted performance ratio >= 26%:")
## [1] "Stores to consider opening based on predicted performance ratio >= 26%:"
print(stores_to_open[, c("STOR", "Performance_Ratio")])
## # A tibble: 2 × 2
##    STOR Performance_Ratio
##   <int>             <dbl>
## 1    57              30.2
## 2    60              36.7


Final Prediction
# Use the complete model (from CroqPainFix) to predict new locations
# Ensure CroqPain_new has the same variable structure and types
# If necessary, convert or create any missing variables

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

# Replace any missing EMPL values with mean (already done earlier)
# Ensure all required variables are present

# Extract the response variable from the final model's formula
response_var <- as.character(formula(model_final)[[2]])

# Get all variables in the formula
all_vars_in_formula <- all.vars(formula(model_final))

# Exclude the response variable to get predictor variables
predictor_vars <- setdiff(all_vars_in_formula, response_var)

# Now, check if all predictor variables are present in CroqPain_new
missing_vars <- setdiff(predictor_vars, names(CroqPain_new))
if (length(missing_vars) > 0) {
  stop("The following variables are missing in CroqPain_new: ", paste(missing_vars, collapse = ", "))
}

# Predict performance for new locations
predictions_new <- predict(model_final, newdata = CroqPain_new)

# Add predictions to CroqPain_new
CroqPain_new$Predicted_EARN <- predictions_new

# Calculate predicted performance ratio
# Assuming performance ratio is EARN / K (investment)
CroqPain_new$Performance_Ratio <- CroqPain_new$Predicted_EARN / CroqPain_new$K * 100  # in percentage

# View predicted performance ratios
print("Predicted Performance Ratios for the new locations:")
## [1] "Predicted Performance Ratios for the new locations:"
print(CroqPain_new[, c("STOR", "Predicted_EARN", "K", "Performance_Ratio")])
##            STOR Predicted_EARN         K Performance_Ratio
## 1        Calais       36.04117  660.1396          5.459629
## 2    Montchanin       60.92370  733.0000          8.311556
## 3      Aubusson       52.57667 1050.3490          5.005638
## 4      Toulouse      267.16341  836.0000         31.957346
## 5         Torcy       43.50262  783.6160          5.551522
## 6  Marseilles-1       92.78076  924.8080         10.032434
## 7  Marseilles-2      -40.40490 1089.5764         -3.708313
## 8      Clermont       81.47771  737.6620         11.045399
## 9   Montpellier      174.29537  584.0000         29.845098
## 10        Dijon      229.11414  681.0000         33.643780
# Recommend which locations to open based on performance ratio target of 26%
locations_to_open <- subset(CroqPain_new, Performance_Ratio >= 26)
print("Recommended locations to open (Performance Ratio >= 26%):")
## [1] "Recommended locations to open (Performance Ratio >= 26%):"
print(locations_to_open[, c("STOR", "Performance_Ratio")])
##           STOR Performance_Ratio
## 4     Toulouse          31.95735
## 9  Montpellier          29.84510
## 10       Dijon          33.64378