Ames Housing Prices Model

library(tidyverse)

library(randomForest)

library(caret)

home_train = read_csv("~/Library/CloudStorage/OneDrive-Personal/Desktop/School/MSDA/Spring 2025/Analytics Applications/Housing Data Project/train.csv")

home_test = read_csv("~/Library/CloudStorage/OneDrive-Personal/Desktop/School/MSDA/Spring 2025/Analytics Applications/Housing Data Project/test.csv")

str(home_train)

anyNA(home_train$SalePrice)

temp_home_test = home_test %>% 
  mutate(SalePrice = NA)
all_homes = rbind(temp_home_test, home_train)

all_homes = all_homes %>% 
  select(-Fireplaces, -MiscFeature)

all_homes = all_homes %>% 
  select(-Id)

identify_categorical = function(data, unique_threshold = 10) {
  is_categorical = sapply(data, function(x) {
    if(is.character(x) || is.factor(x) || is.logical(x)) {
      return(TRUE)
    }
    
    if(is.numeric(x) && length(unique(x)) <= unique_threshold) {
      return(TRUE)
    }
    
    return(FALSE)
  })
  
  return(names(data)[is_categorical])
}

#Identify categorical variables
categorical_vars = identify_categorical(all_homes)
cat("Identified categorical variables:\n")

print(categorical_vars)

#Plot
plot_categorical = function(data, var_names) {
  plots = list()
  
  for(var in var_names) {
    if(!is.factor(data[[var]])) {
      data[[var]] = as.factor(data[[var]])
    }
    
    var_counts = data %>%
      count(!!sym(var)) %>%
      mutate(percentage = n / sum(n) * 100) %>%
      arrange(desc(n))
    
    p = ggplot(var_counts, aes(x = reorder(!!sym(var), -n), y = n)) +
      geom_bar(stat = "identity", fill = "steelblue", width = 0.7) +
      geom_text(aes(label = paste0(n, " (", round(percentage, 1), "%)")), 
                vjust = -0.5, size = 3) +
      labs(title = paste("Distribution of", var),
           subtitle = paste(length(unique(data[[var]])), "unique values"),
           x = var,
           y = "Count") +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1),
            plot.title = element_text(face = "bold"),
            panel.grid.minor = element_blank())
    
    if(length(unique(data[[var]])) > 15) {
      p = p + scale_x_discrete(guide = guide_axis(n.dodge = 2))
      cat(paste0("Note: Variable '", var, "' has many categories (", 
                 length(unique(data[[var]])), "). Consider grouping some categories.\n"))
    }
    
    plots[[var]] = p
    
    print(p)
  }
  
  return(plots)
}

cat_plots = plot_categorical(all_homes, categorical_vars)

colSums(is.na(all_homes))

all_homes = all_homes %>% 
  select(-Fence, -Alley, -PoolQC, -PoolArea, -Street, -Utilities, -LandSlope, -Condition2, -RoofMatl, -Heating, -KitchenAbvGr)

any(duplicated(all_homes) == TRUE)

# Function to identify continuous variables
identify_continuous = function(data, unique_threshold = 10) {
  is_continuous = sapply(data, function(x) {
    # Check if numeric and has more unique values than our threshold
    if(is.numeric(x) && length(unique(x)) > unique_threshold) {
      return(TRUE)
    }
    return(FALSE)
  })
  
  return(names(data)[is_continuous])
}

# Identify continuous variables
continuous_vars = identify_continuous(all_homes)
cat("Identified continuous variables:\n")

print(continuous_vars)

# Create histograms for each continuous variable
plot_histograms = function(data, var_names) {
  plots = list()
  
  for(var in var_names) {
    # Generate basic histogram
    p = ggplot(data, aes(x = !!sym(var))) +
      geom_histogram(fill = "steelblue", color = "white", bins = 30) +
      labs(title = paste("Distribution of", var),
           subtitle = paste("n =", sum(!is.na(data[[var]])), 
                           "| NA =", sum(is.na(data[[var]]))),
           x = var,
           y = "Count") +
      theme_minimal() +
      theme(plot.title = element_text(face = "bold"),
            panel.grid.minor = element_blank())
    
    # Add a vertical line for the mean
    p = p + geom_vline(aes(xintercept = mean(!!sym(var), na.rm = TRUE)),
                       color = "darkred", linetype = "dashed", linewidth = 1)
    
    # Save plot to list
    plots[[var]] = p
    
    # Display plot
    print(p)
  }
  
  return(plots)
}

# Generate and display all histograms
hist_plots = plot_histograms(all_homes, continuous_vars)

# If you want to save all plots to PDF
save_histograms = function(plot_list, filename = "continuous_histograms.pdf") {
  pdf(filename, width = 10, height = 7)
  for(p in plot_list) {
    print(p)
  }
  dev.off()
  cat(paste0("All histograms saved to ", filename, "\n"))
}

save_histograms(hist_plots)

which(all_homes$GarageYrBlt == max(range(all_homes$GarageYrBlt, na.rm = T)))

max(range(all_homes$GarageYrBlt, na.rm = T))

all_homes$GarageYrBlt[1133] = 2007
all_homes$GarageYrBlt[1133]

outliers = c()

outliers = append(outliers, which(all_homes$MiscVal > 100))
outliers = append(outliers, which(all_homes$LotArea > 50000))
outliers = append(outliers, which(all_homes$LotFrontage == max(range(all_homes$LotFrontage, na.rm = T))))
outliers = append(outliers, which(all_homes$MasVnrArea > 1000))
outliers = append(outliers, which(all_homes$BsmtFinSF1 == max(range(all_homes$BsmtFinSF1, na.rm = T))))
outliers = append(outliers, which(all_homes$BsmtFinSF2 == max(range(all_homes$BsmtFinSF2, na.rm = T))))
outliers = append(outliers, which(all_homes$TotalBsmtSF == max(range(all_homes$TotalBsmtSF, na.rm = T))))
outliers = append(outliers, which(all_homes$`1stFlrSF` == max(range(all_homes$`1stFlrSF`, na.rm = T))))
outliers = append(outliers, which(all_homes$`2ndFlrSF` == max(range(all_homes$`2ndFlrSF`, na.rm = T))))
outliers = append(outliers, which(all_homes$LowQualFinSF > 50))
outliers = append(outliers, which(all_homes$GrLivArea == max(range(all_homes$GrLivArea, na.rm = T))))
outliers = append(outliers, which(all_homes$WoodDeckSF == max(range(all_homes$WoodDeckSF, na.rm = T))))
outliers = append(outliers, which(all_homes$OpenPorchSF == max(range(all_homes$OpenPorchSF, na.rm = T))))
outliers = append(outliers, which(all_homes$EnclosedPorch > 400))
outliers = append(outliers, which(all_homes$`3SsnPorch` > 100))
outliers = append(outliers, which(all_homes$ScreenPorch == max(range(all_homes$ScreenPorch, na.rm = T))))
outliers = unique(outliers)
all_homes = all_homes[-outliers, ]

hist_plots = plot_histograms(all_homes, continuous_vars)

all_homes = all_homes %>% 
  select(-MiscVal, -ScreenPorch, -`3SsnPorch`, -EnclosedPorch, -LowQualFinSF, -BsmtFinSF2)

garage_bsmt_cols = grep("Garage|Bsmt", names(all_homes), value = TRUE)

all_homes = all_homes %>%
  mutate(across(all_of(garage_bsmt_cols), 
                ~if(is.numeric(.)) replace_na(., 0) else 
                  if(is.character(.)) replace_na(., "0") else 
                    if(is.factor(.)) {
                      f = .
                      if(!("0" %in% levels(f))) {
                        levels(f) = c(levels(f), "0")
                      }
                      f[is.na(f)] = "0"
                      f
                    } else .))

na_counts = colSums(is.na(all_homes[, garage_bsmt_cols, drop = FALSE]))
print(na_counts)

colSums(is.na(all_homes))

all_homes$FireplaceQu = replace_na(all_homes$FireplaceQu, "0")

get_mode = function(x) {
  # Remove NAs for calculation
  x = x[!is.na(x)]
  
  # If empty vector, return NA
  if(length(x) == 0) return(NA)
  
  # Find most frequent value
  unique_values = unique(x)
  unique_values[which.max(tabulate(match(x, unique_values)))]
}

# Function to replace NAs in the dataframe
replace_nas = function(data) {
  result = data
  
  for(col in names(data)) {
    # Skip if no NAs
    if(!any(is.na(data[[col]]))) next
    
    # For numeric columns, replace with mean
    if(is.numeric(data[[col]])) {
      mean_val = mean(data[[col]], na.rm = TRUE)
      result[[col]][is.na(result[[col]])] = mean_val
      cat(paste("Replaced NAs in numeric column", col, "with mean:", mean_val, "\n"))
    }
    # For character columns, replace with mode
    else if(is.character(data[[col]])) {
      mode_val = get_mode(data[[col]])
      result[[col]][is.na(result[[col]])] = mode_val
      cat(paste("Replaced NAs in character column", col, "with mode:", mode_val, "\n"))
    }
    # For factor columns, also replace with mode
    else if(is.factor(data[[col]])) {
      mode_val = get_mode(as.character(data[[col]]))
      # Make sure mode is in the levels
      if(!mode_val %in% levels(result[[col]])) {
        levels(result[[col]]) = c(levels(result[[col]]), mode_val)
      }
      result[[col]][is.na(result[[col]])] = mode_val
      cat(paste("Replaced NAs in factor column", col, "with mode:", mode_val, "\n"))
    }
  }
  
  return(result)
}

all_homes_clean = all_homes %>% 
  select(-SalePrice) %>% 
  replace_nas()

na_count = sapply(all_homes_clean, function(x) sum(is.na(x)))
if(any(na_count > 0)) {
  cat("\nRemaining NAs in columns:\n")
  print(na_count[na_count > 0])
} else {
  cat("\nSuccess! All NAs have been replaced.\n")
}

all_homes_clean$SalePrice = all_homes$SalePrice
colSums(is.na(all_homes_clean))

all_homes_clean$perc1stSF = all_homes_clean$`1stFlrSF`/all_homes_clean$GrLivArea
all_homes_clean$perc2ndSF = all_homes_clean$`2ndFlrSF`/all_homes_clean$GrLivArea
all_homes_clean = all_homes_clean %>% 
  select(-`1stFlrSF`, -`2ndFlrSF`)

training_index = is.na(all_homes_clean$SalePrice)
train = all_homes_clean[training_index == F, ]
test = all_homes_clean[training_index == T, ]
test = test %>% 
  select(-SalePrice)

dim(train)

dim(test)

train %>% 
  group_by(FireplaceQu) %>% 
  summarise(Count = n())

fire = train %>%
  mutate(FireplaceStatus = ifelse(FireplaceQu == "0", "No Fireplace", "Has Fireplace"))

fire %>% 
  ggplot(aes(x = FireplaceStatus, y = SalePrice/1000, fill = FireplaceStatus, colour = FireplaceStatus)) +
  geom_boxplot(alpha = 0.5) +
  scale_y_continuous(labels = scales::label_dollar()) +
  xlab("Fireplace Status") +
  ylab("Sale Price in 1000s") +
  ggtitle("Are Fireplaces Adding Value to Homes?") +
  theme_minimal() +
  theme(legend.position = "none")

train %>% 
  ggplot(aes(x = factor(FireplaceQu, levels = c("Ex", "Gd", "TA", "Fa", "Po", "0"), ordered = TRUE), y = SalePrice/1000, fill = FireplaceQu, colour = FireplaceQu)) +
  geom_boxplot(alpha = 0.5) +
  scale_y_continuous(labels = scales::label_dollar()) +
  xlab("Fireplace Quality") +
  ylab("Sale Price in 1000s") +
  ggtitle("Are Quality Fireplaces Adding Value to Homes?") +
  theme_minimal() +
  theme(legend.position = "none") 

all_homes_clean %>% 
  group_by(Functional) %>% 
  summarise(Count = n())

all_homes_clean %>% 
  ggplot(aes(x = factor(Functional, levels = c("Typ", "Min1", "Min2", "Mod", "Maj1", "Maj2", "Sev"), ordered = TRUE), y = YearBuilt, fill = Functional, colour = Functional)) +
  geom_boxplot(alpha = 0.5) +
  xlab("Functionality") +
  ylab("Year Built") +
  ggtitle("Are Older Homes Considered Less Functional?") +
  theme_minimal() +
  theme(legend.position = "none") 

all_homes_clean %>% 
  group_by(HalfBath) %>% 
  summarise(Count = n())

stories = all_homes_clean %>%
  mutate(numStories = ifelse(perc2ndSF == "0", "One Story", "Two Story"))

stories %>% 
  ggplot(aes(fill = numStories, x = HalfBath)) +
    geom_bar(position = "dodge", alpha = 0.5) +
    xlab("Half Baths") +
    ylab("Count") +
    ggtitle("Are 2-Story Homes More Likely to have a Half Bath?") +
    labs(fill = "Number of Stories") +
    theme_minimal() 

set.seed(23)

trainIndex = createDataPartition(train$SalePrice, p = 0.7, list = FALSE)
train_set = train[trainIndex, ]
validation_set = train[-trainIndex, ]

cat("Training set size:", nrow(train_set), "\n")

cat("Validation set size:", nrow(validation_set), "\n")

set.seed(23)

rf_model = randomForest(
  x = train_set[, setdiff(names(train_set), "SalePrice")],
  y = train_set$SalePrice,
  importance = TRUE,
  ntree = 500
)

importance_df = as.data.frame(randomForest::importance(rf_model))
importance_df$Variable = rownames(importance_df)

sorted_importance = importance_df[order(importance_df$`%IncMSE`, decreasing = TRUE), ]

top_10_vars = head(sorted_importance$Variable, 10)

cat("Top 10 most important variables:\n")

print(top_10_vars)

sorted_importance[1:10, ] %>%
  ggplot(aes(reorder(x = Variable, `%IncMSE`), y = `%IncMSE`)) +
  geom_bar(stat = "identity", fill = "orange", color = "black", width = 0.5) +
  coord_flip() +
  labs(x = "Variables", y = "Variable importance") +
  theme_minimal()

train_set_top10 = train_set[, c(top_10_vars, "SalePrice")]
validation_set_top10 = validation_set[, c(top_10_vars, "SalePrice")]

cat("Dimensions of reduced training set:", dim(train_set_top10), "\n")

cat("Dimensions of reduced validation set:", dim(validation_set_top10), "\n")

library(xgboost)

library(e1071)    
library(BART)       

library(car)       

library(corrplot) 

library(ranger)

set.seed(23)

ctrl = trainControl(
  method = "cv",           
  number = 5,              
  verboseIter = TRUE,      
  savePredictions = "final" 
)


build_and_evaluate = function(model_type, train_data, test_data, tune_grid = NULL) {
  formula = as.formula("SalePrice ~ .")
  
  # Train model
  if (!is.null(tune_grid)) {
    model = train(
      formula,
      data = train_data,
      method = model_type,
      trControl = ctrl,
      metric = "RMSE",
      tuneGrid = tune_grid
    )
  } else {
    model = train(
      formula,
      data = train_data,
      method = model_type,
      trControl = ctrl,
      metric = "RMSE"
    )
  }
  
  # Make predictions
  predictions = predict(model, newdata = test_data)
  
  # Calculate RMSE
  rmse = sqrt(mean((predictions - test_data$SalePrice)^2))
  
  # Calculate R-squared
  rsq = 1 - sum((test_data$SalePrice - predictions)^2) / 
    sum((test_data$SalePrice - mean(test_data$SalePrice))^2)
  
  # Return results
  return(list(
    model = model,
    predictions = predictions,
    rmse = rmse,
    rsq = rsq
  ))
}

# Random Forest
rf_grid = expand.grid(
  mtry = c(3, 4, 5, 6)
)

# SVM
svm_grid = expand.grid(
  sigma = c(0.01, 0.1, 1),
  C = c(1, 10, 100)
)

# XGBoost
xgb_grid = expand.grid(
  nrounds = c(100, 200),
  max_depth = c(3, 6),
  eta = c(0.1, 0.3),
  gamma = 0,
  colsample_bytree = 0.8,
  min_child_weight = 1,
  subsample = 1
)

linear_train_top10 = train_set_top10 %>% 
  mutate(ExterQual = as.factor(ExterQual))
linear_val_top10 = validation_set_top10 %>% 
  mutate(ExterQual = as.factor(ExterQual))
correlation_matrix = linear_train_top10 %>% 
  select(-ExterQual) %>% 
  cor()
print("Correlation matrix of predictors:")

print(round(correlation_matrix, 2))

corrplot(correlation_matrix, method = "color", type = "upper", 
         tl.col = "black", tl.srt = 45)

initial_lm = lm(SalePrice ~ ., data = linear_train_top10)
vif_values = vif(initial_lm)
print("Variance Inflation Factors (VIF):")

print(vif_values)

models = list()
models$lm = build_and_evaluate("lm", linear_train_top10, linear_val_top10)

set.seed(23)

# Random Forest
models$rf = build_and_evaluate("rf", train_set_top10, validation_set_top10, rf_grid)

# SVM
models$svm = build_and_evaluate("svmRadial", train_set_top10, validation_set_top10, svm_grid)

# XGBoost
models$xgb = build_and_evaluate("xgbTree", train_set_top10, validation_set_top10, xgb_grid)

# Compare results
results = data.frame(
  Model = character(),
  RMSE = numeric(),
  Rsquared = numeric(),
  stringsAsFactors = FALSE
)

for (name in names(models)) {
  if (!is.null(models[[name]])) {
    results = rbind(results, data.frame(
      Model = name,
      RMSE = models[[name]]$rmse,
      Rsquared = models[[name]]$rsq
    ))
  }
}

print(results[order(results$RMSE), ])

library(glmnet)     # For ridge and lasso regression
library(rpart)      # For decision trees
library(gbm)        # For gradient boosting
library(earth)      # For MARS (Multivariate Adaptive Regression Splines)
library(Cubist)     # For Cubist models
library(kernlab)    # For kernlab SVM implementation

set.seed(23)

# 1. Ridge Regression
ridge_grid = expand.grid(alpha = 0, 
                         lambda = 10^seq(-3, 3, length = 10))
models$ridge = build_and_evaluate("glmnet", train_set_top10, validation_set_top10, ridge_grid)

# 2. Lasso Regression
lasso_grid = expand.grid(alpha = 1, 
                         lambda = 10^seq(-3, 3, length = 10))
models$lasso = build_and_evaluate("glmnet", train_set_top10, validation_set_top10, lasso_grid)

# 3. Elastic Net (combination of Ridge and Lasso)
elasticnet_grid = expand.grid(alpha = seq(0, 1, length = 5), 
                              lambda = 10^seq(-3, 3, length = 6))
models$elasticnet = build_and_evaluate("glmnet", train_set_top10, validation_set_top10, elasticnet_grid)

# 4. Decision Tree with tuning
rpart_grid = expand.grid(cp = seq(0.001, 0.1, length = 10))
models$rpart = build_and_evaluate("rpart", train_set_top10, validation_set_top10, rpart_grid)

# 5. Gradient Boosting Machine (GBM)
gbm_grid = expand.grid(
  n.trees = c(100, 500),
  interaction.depth = c(3, 5, 7),
  shrinkage = c(0.01, 0.1),
  n.minobsinnode = 10
)
models$gbm = build_and_evaluate("gbm", train_set_top10, validation_set_top10, gbm_grid)

# 6. MARS (Multivariate Adaptive Regression Splines)
mars_grid = expand.grid(
  degree = 1:3,
  nprune = c(5, 10, 15, 20)
)
models$mars = build_and_evaluate("earth", train_set_top10, validation_set_top10, mars_grid)

# 7. Cubist model (rule-based model)
cubist_grid = expand.grid(
  committees = c(1, 10, 50, 100),
  neighbors = c(0, 5, 9)
)
models$cubist = build_and_evaluate("cubist", train_set_top10, validation_set_top10, cubist_grid)

# 8. Gaussian Process
models$gausspr = build_and_evaluate("gaussprRadial", train_set_top10, validation_set_top10)

# 9. k-Nearest Neighbors
knn_grid = expand.grid(k = c(3, 5, 7, 9))
models$knn = build_and_evaluate("knn", train_set_top10, validation_set_top10, knn_grid)

# Update results table with all models
results = data.frame(
  Model = character(),
  RMSE = numeric(),
  Rsquared = numeric(),
  stringsAsFactors = FALSE
)

for (name in names(models)) {
  if (!is.null(models[[name]])) {
    results = rbind(results, data.frame(
      Model = name,
      RMSE = models[[name]]$rmse,
      Rsquared = models[[name]]$rsq
    ))
  }
}

print(results[order(results$RMSE), ])

ggplot(results, aes(x = reorder(Model, -RMSE), y = RMSE)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = round(RMSE, 0)), hjust = -0.1) +
  coord_flip(ylim = c(20000,45000)) +
  labs(title = "Model Comparison by RMSE",
       x = "Model",
       y = "RMSE (lower is better)") +
  theme_minimal()

ggplot(results, aes(x = reorder(Model, Rsquared), y = Rsquared)) +
  geom_bar(stat = "identity", fill = "darkgreen") +
  geom_text(aes(label = round(Rsquared, 3)), hjust = -0.1) +
  coord_flip(ylim = c(0.65, 1.0)) +
  labs(title = "Model Comparison by R-squared",
       x = "Model",
       y = "R-squared (higher is better)") +
  theme_minimal()

top3_models = results$Model[order(results$RMSE)][1:3]
cat("Top 3 models:", paste(top3_models, collapse=", "), "\n")

set.seed(23)

# Create average predictions
ensemble_preds = numeric(nrow(validation_set_top10))
for (model_name in top3_models) {
  preds = models[[model_name]]$predictions
  if (length(preds) != nrow(validation_set_top10)) {
    stop(paste("Predictions from", model_name, "are not the right length"))
  }
  print(paste("Adding predictions from model:", model_name))
  ensemble_preds = ensemble_preds + preds
}

ensemble_preds = ensemble_preds / length(top3_models)


# Calculate performance metrics for ensemble
ensemble_rmse = sqrt(mean((ensemble_preds - validation_set_top10$SalePrice)^2))
ensemble_rsq = 1 - sum((validation_set_top10$SalePrice - ensemble_preds)^2) / 
  sum((validation_set_top10$SalePrice - mean(validation_set_top10$SalePrice))^2)

# Add ensemble to results
results = rbind(results, data.frame(
  Model = "ensemble_top3",
  RMSE = ensemble_rmse,
  Rsquared = ensemble_rsq
))

print(results[order(results$RMSE), ])

plot_data = data.frame(
  Actual = validation_set_top10$SalePrice,
  Predicted = ensemble_preds
)

ggplot(plot_data, aes(x = Actual, y = Predicted)) +
  geom_point(color = "steelblue", alpha = 0.6) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Actual vs. Predicted Sale Prices (Ensemble)",
    x = "Actual Sale Price",
    y = "Predicted Sale Price"
  ) +
  theme_minimal()

top3_models

set.seed(23)

predictions1 = predict(models$svm$model, newdata = test)
predictions2 = predict(models$xgb$model, newdata = test)
predictions3 = predict(models$gbm$model, newdata = test)

avg_test_preds = (predictions1 + predictions2 + predictions3)/3
test$SalePricePred = round(avg_test_preds,-2)

head(test[,57:61])