1. Introduction

The New York City property market is large and active. Every year, many buildings, houses, and apartment units are sold across different boroughs and neighborhoods. Property sale price may be affected by several factors, such as location, building class, land size, gross square feet, number of units, year built, and sale date.

In this project, we use the NYC Property Sales dataset to study property transaction patterns in New York City. The dataset records property sales over a 12-month period and includes information about location, property type, property size, number of units, sale price, and sale date.

The purpose of this project is to use R programming to clean, explore, visualize, and model the dataset. The project focuses on two main tasks: property sale price prediction and transaction classification. The price prediction task is treated as a regression problem, while the transaction classification task is treated as a classification problem.

2. Project Objectives and Research Questions

The project has three main objectives:

To explore the main property features that are related to sale price in New York City.
To predict property sale price using selected property features through a regression approach.
To classify property transactions into meaningful groups, such as High Value and Regular Value transactions.

This project aims to answer the following research questions:

Research Question 1: What are the main property features related to sale price in New York City?

This question helps us understand whether variables such as borough, neighborhood, building class, land square feet, gross square feet, total units, and year built are related to sale price.

Research Question 2: Can we predict property sale price using property features?

This is a regression problem. The target variable is SALE PRICE.

Research Question 3: Can we classify whether a property transaction is a High Value transaction?

This is a classification problem. A new target variable can be created based on sale price. For example, transactions above a selected price threshold can be classified as High Value, while the remaining transactions can be classified as Regular Value.

3. Dataset Description

Dataset Source, Period, and Purpose

The dataset used in this project is the NYC Property Sales dataset from Kaggle.

Dataset link: https://www.kaggle.com/datasets/new-york-city/nyc-property-sales

The original data comes from the New York City Department of Finance Rolling Sales Data. The dataset records property sales in New York City over a 12-month period. Based on the SALE DATE variable, the dataset mainly covers property transactions from 2016 to 2017.

The main purpose of using this dataset is to understand NYC property sales patterns and support two modelling tasks: sale price prediction and transaction value classification.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readr)
library(knitr)

nyc_sales <- read_csv("nyc-rolling-sales.csv", show_col_types = FALSE)

## New names:
## • `` -> `...1`

Dataset Dimension

The raw dataset contains 84,548 rows and 22 columns.

raw_rows <- nrow(nyc_sales)
raw_columns <- ncol(nyc_sales)
raw_dimension <- raw_rows * raw_columns

dimension_table <- data.frame(
  Item = c("Number of rows", "Number of columns", "Total dimension"),
  Value = c(raw_rows, raw_columns, raw_dimension)
)

kable(dimension_table, row.names = FALSE)

Item	Value
Number of rows	84548
Number of columns	22
Total dimension	1860056

The raw dataset dimension is:

84,548 rows × 22 columns = 1,860,056

This is much larger than 100,000, so the dataset satisfies the project requirement.

The first column is an unnecessary index column from the original file. After removing this column, the dataset has 21 named data columns for further checking and analysis.

nyc_sales_start <- nyc_sales %>%
  select(-1)

useful_rows <- nrow(nyc_sales_start)
useful_columns <- ncol(nyc_sales_start)
useful_dimension <- useful_rows * useful_columns

useful_dimension_table <- data.frame(
  Item = c("Number of rows", "Number of data columns after removing index", "Dataset dimension after removing index"),
  Value = c(useful_rows, useful_columns, useful_dimension)
)

kable(useful_dimension_table, row.names = FALSE)

Item	Value
Number of rows	84548
Number of data columns after removing index	21
Dataset dimension after removing index	1775508

After removing the unnecessary index column, the dataset dimension is:

84,548 rows × 21 columns = 1,775,508

This still satisfies the project requirement.

Key Variables and Data Structure

Each row in the dataset represents one property sale transaction in New York City. Each column describes one feature of the property or the sale transaction.

variable_structure <- data.frame(
  Variable = names(nyc_sales_start),
  Data_Type = vapply(
    nyc_sales_start,
    function(x) paste(class(x), collapse = ", "),
    character(1)
  )
)

kable(variable_structure, row.names = FALSE)

Variable	Data_Type
BOROUGH	numeric
NEIGHBORHOOD	character
BUILDING CLASS CATEGORY	character
TAX CLASS AT PRESENT	character
BLOCK	numeric
LOT	numeric
EASE-MENT	logical
BUILDING CLASS AT PRESENT	character
ADDRESS	character
APARTMENT NUMBER	character
ZIP CODE	numeric
RESIDENTIAL UNITS	numeric
COMMERCIAL UNITS	numeric
TOTAL UNITS	numeric
LAND SQUARE FEET	character
GROSS SQUARE FEET	character
YEAR BUILT	numeric
TAX CLASS AT TIME OF SALE	numeric
BUILDING CLASS AT TIME OF SALE	character
SALE PRICE	character
SALE DATE	POSIXct, POSIXt

The dataset includes location variables, property type variables, unit variables, size variables, time variables, and the target variable for price prediction.

Variable Group	Variables	Description
Location variables	BOROUGH, NEIGHBORHOOD, ZIP CODE, ADDRESS	Show where the property is located
Property type variables	BUILDING CLASS CATEGORY, TAX CLASS AT PRESENT, BUILDING CLASS AT PRESENT	Describe property type and tax class
Unit variables	RESIDENTIAL UNITS, COMMERCIAL UNITS, TOTAL UNITS	Show the number of residential, commercial, and total units
Size variables	LAND SQUARE FEET, GROSS SQUARE FEET	Show land size and building size
Time variables	YEAR BUILT, SALE DATE	Show building year and transaction date
Target variable	SALE PRICE	Shows the final sale price of the property

The most important variables for this project are shown below:

Variable	Description	Use in This Project
BOROUGH	The borough where the property is located	Location analysis and modelling feature
NEIGHBORHOOD	More detailed location information	Price comparison by area
BUILDING CLASS CATEGORY	General property type	Property type analysis
RESIDENTIAL UNITS	Number of residential units	Modelling feature
COMMERCIAL UNITS	Number of commercial units	Modelling feature
TOTAL UNITS	Total number of units	Modelling feature
LAND SQUARE FEET	Land area of the property	Price prediction feature
GROSS SQUARE FEET	Total building area	Price prediction feature
YEAR BUILT	Year when the property was built	Property age analysis
SALE DATE	Date of property sale	Time trend analysis
SALE PRICE	Final sale price	Target variable for regression

Initial Dataset Summary

Overall, the NYC Property Sales dataset is suitable for this project because it is large enough and contains both numerical and categorical variables. It can support data cleaning, exploratory data analysis, visualization, regression modelling, and classification modelling.

For Research Question 1, the dataset can be used to explore which property features are related to sale price. For Research Question 2, SALE PRICE can be used as the target variable for regression. For Research Question 3, a new classification target can be created to classify property transactions into High Value and Regular Value groups based on sale price.

However, the raw dataset has some data quality issues. For example, SALE PRICE, LAND SQUARE FEET, and GROSS SQUARE FEET may need type conversion before modelling. Some missing or invalid values may also appear in the dataset. In addition, the first column is only an index column, and some columns need to be checked carefully during the data cleaning stage. These issues will be handled in the data cleaning and preprocessing section before analysis and modelling.

4. Data Cleaning and Preprocessing

The raw data required cleaning before exploratory analysis and modelling. The main issues included an unnecessary index column, a non-informative EASE-MENT field, numeric values stored as text, missing values, duplicate records, invalid or non-market sale prices, unrealistic property sizes, and invalid building years.

Loading the Raw Data

raw_file_candidates <- c(
  "nyc-rolling-sales.csv",
  "data/nyc-rolling-sales.csv",
  "../data/nyc-rolling-sales.csv",
  "work/nyc_property_sales_raw/nyc-rolling-sales.csv"
)

raw_file <- raw_file_candidates[file.exists(raw_file_candidates)][1]

if (is.na(raw_file)) {
  stop("Raw file not found. Please place Kaggle's nyc-rolling-sales.csv in the same folder or in a data/ folder.")
}

raw_data <- read.csv(
  raw_file,
  stringsAsFactors = FALSE,
  check.names = FALSE,
  na.strings = c("", " ", "-", " -  ", "NA", "N/A")
)

raw_rows <- nrow(raw_data)
raw_cols <- ncol(raw_data)
cat("Raw dataset dimension:", raw_rows, "rows x", raw_cols, "columns\n")

## Raw dataset dimension: 84548 rows x 22 columns

Cleaning Column Names and Values

clean_column_name <- function(x) {
  x <- tolower(trimws(x))
  x <- gsub("[^a-z0-9]+", "_", x)
  x <- gsub("^_+|_+$", "", x)
  x
}

clean_text <- function(x) {
  x <- trimws(as.character(x))
  x[x %in% c("", "-", "NA", "N/A", "NULL", "null")] <- NA
  x
}

to_number <- function(x) {
  as.numeric(gsub(",", "", clean_text(x)))
}

data <- raw_data
names(data) <- clean_column_name(names(data))

# The first Kaggle column is only a row index, and EASE-MENT contains no useful information.
# R may read the blank Kaggle index header as "" in some environments.
drop_cols <- intersect(c("", "x", "unnamed", "unnamed_0", "ease_ment"), names(data))
data <- data[, !(names(data) %in% drop_cols)]

character_cols <- names(data)[sapply(data, is.character)]
for (col in character_cols) {
  data[[col]] <- clean_text(data[[col]])
}

cat("Columns after dropping index/non-informative fields:", ncol(data), "\n")

## Columns after dropping index/non-informative fields: 20

Data Type Conversion

numeric_from_text <- intersect(
  c("sale_price", "land_square_feet", "gross_square_feet"),
  names(data)
)

for (col in numeric_from_text) {
  data[[col]] <- to_number(data[[col]])
}

integer_cols <- intersect(
  c(
    "borough", "block", "lot", "zip_code", "residential_units",
    "commercial_units", "total_units", "year_built",
    "tax_class_at_time_of_sale"
  ),
  names(data)
)

for (col in integer_cols) {
  data[[col]] <- as.integer(data[[col]])
}

data$sale_date <- as.Date(substr(data$sale_date, 1, 10))
data$sale_year <- as.integer(format(data$sale_date, "%Y"))
data$sale_month <- as.integer(format(data$sale_date, "%m"))

borough_lookup <- c(
  "1" = "Manhattan",
  "2" = "Bronx",
  "3" = "Brooklyn",
  "4" = "Queens",
  "5" = "Staten Island"
)

data$borough_name <- borough_lookup[as.character(data$borough)]
data$borough_name[is.na(data$borough_name)] <- "Unknown"

if ("zip_code" %in% names(data)) {
  data$zip_code[data$zip_code == 0] <- NA
}

Missing Values and Invalid Records

cleaning_log <- data.frame(
  step = "Raw data",
  rows = raw_rows,
  columns = raw_cols
)

add_log <- function(log_data, step_name, current_data) {
  rbind(
    log_data,
    data.frame(step = step_name, rows = nrow(current_data), columns = ncol(current_data))
  )
}

# Fill categorical missing values with "Unknown".
categorical_cols <- intersect(
  c(
    "neighborhood", "building_class_category", "tax_class_at_present",
    "building_class_at_present", "address", "apartment_number",
    "building_class_at_time_of_sale", "borough_name"
  ),
  names(data)
)

for (col in categorical_cols) {
  data[[col]][is.na(data[[col]])] <- "Unknown"
}

# Median imputation for optional numeric land size.
if ("land_square_feet" %in% names(data)) {
  land_median <- median(data$land_square_feet, na.rm = TRUE)
  data$land_square_feet[is.na(data$land_square_feet)] <- land_median
}

data <- unique(data)
cleaning_log <- add_log(cleaning_log, "Remove duplicate rows", data)

# Required fields for price prediction and feature engineering.
data <- data[!is.na(data$sale_date), ]
cleaning_log <- add_log(cleaning_log, "Keep records with valid sale date", data)

data <- data[!is.na(data$sale_price) & data$sale_price >= 10000, ]
cleaning_log <- add_log(cleaning_log, "Keep market transactions with sale_price >= 10000", data)

data <- data[!is.na(data$gross_square_feet) & data$gross_square_feet >= 200, ]
cleaning_log <- add_log(cleaning_log, "Keep records with gross_square_feet >= 200", data)

data <- data[
  !is.na(data$year_built) &
    data$year_built >= 1800 &
    !is.na(data$sale_year) &
    data$year_built <= data$sale_year,
]
cleaning_log <- add_log(cleaning_log, "Keep valid building years", data)

data <- data[
  !is.na(data$residential_units) &
    !is.na(data$commercial_units) &
    !is.na(data$total_units) &
    data$residential_units >= 0 &
    data$commercial_units >= 0 &
    data$total_units >= 0,
]
cleaning_log <- add_log(cleaning_log, "Keep valid unit counts", data)

Cleaned Dataset Summary

cleaning_log

##                                                     step  rows columns
## 1                                               Raw data 84548      22
## 2                                  Remove duplicate rows 83783      23
## 3                      Keep records with valid sale date 83783      23
## 4      Keep market transactions with sale_price >= 10000 58333      23
## 5             Keep records with gross_square_feet >= 200 28333      23
## 6                              Keep valid building years 28323      23
## 7                                 Keep valid unit counts 28323      23
## 8                Keep price_per_sqft between 10 and 5000 28123      25
## 9         Remove sale_price outside 1st-99th percentiles 27562      25
## 10 Remove gross_square_feet outside 1st-99th percentiles 27025      25
## 11        Final cleaned dataset with engineered features 27025      27

summary(cleaned_data[, c("sale_price", "gross_square_feet", "property_age", "price_per_sqft")])

##    sale_price       gross_square_feet  property_age    price_per_sqft   
##  Min.   :   90000   Min.   :  720     Min.   :  0.00   Min.   :  11.89  
##  1st Qu.:  446160   1st Qu.: 1360     1st Qu.: 57.00   1st Qu.: 251.03  
##  Median :  635000   Median : 1848     Median : 86.00   Median : 352.68  
##  Mean   :  963664   Mean   : 2346     Mean   : 75.88   Mean   : 409.63  
##  3rd Qu.:  950000   3rd Qu.: 2580     3rd Qu.: 97.00   3rd Qu.: 493.01  
##  Max.   :16232000   Max.   :21750     Max.   :217.00   Max.   :4981.88

table(cleaned_data$borough_name)

## 
##         Bronx      Brooklyn     Manhattan        Queens Staten Island 
##          3200          7988           640         10453          4744

table(cleaned_data$high_value_label)

## 
##    High Value Regular Value 
##          6888         20137

Exporting the Cleaned Dataset

output_dir <- "outputs"
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

write.csv(
  cleaned_data,
  file.path(output_dir, "nyc_property_sales_cleaned.csv"),
  row.names = FALSE
)

write.csv(
  cleaning_log,
  file.path(output_dir, "nyc_property_sales_cleaning_log.csv"),
  row.names = FALSE
)

cat("Cleaned data exported to:", file.path(output_dir, "nyc_property_sales_cleaned.csv"), "\n")

## Cleaned data exported to: outputs/nyc_property_sales_cleaned.csv

cat("Cleaning log exported to:", file.path(output_dir, "nyc_property_sales_cleaning_log.csv"), "\n")

## Cleaning log exported to: outputs/nyc_property_sales_cleaning_log.csv

5. Feature Engineering

data$property_age <- data$sale_year - data$year_built
data$price_per_sqft <- data$sale_price / data$gross_square_feet

# Remove impossible or extreme price-per-square-foot values.
data <- data[
  !is.na(data$price_per_sqft) &
    data$price_per_sqft >= 10 &
    data$price_per_sqft <= 5000,
]
cleaning_log <- add_log(cleaning_log, "Keep price_per_sqft between 10 and 5000", data)

keep_middle_percentile <- function(df, col, lower = 0.01, upper = 0.99) {
  limits <- quantile(df[[col]], probs = c(lower, upper), na.rm = TRUE)
  df[df[[col]] >= limits[1] & df[[col]] <= limits[2], ]
}

# Trim only the most extreme 1% tails for model stability.
data <- keep_middle_percentile(data, "sale_price")
cleaning_log <- add_log(cleaning_log, "Remove sale_price outside 1st-99th percentiles", data)

data <- keep_middle_percentile(data, "gross_square_feet")
cleaning_log <- add_log(cleaning_log, "Remove gross_square_feet outside 1st-99th percentiles", data)

high_value_threshold <- as.numeric(quantile(data$sale_price, 0.75, na.rm = TRUE))
data$high_value <- ifelse(data$sale_price >= high_value_threshold, 1, 0)
data$high_value_label <- ifelse(data$high_value == 1, "High Value", "Regular Value")

factor_cols <- intersect(
  c(
    "borough_name", "neighborhood", "building_class_category",
    "tax_class_at_present", "building_class_at_present",
    "building_class_at_time_of_sale", "high_value_label"
  ),
  names(data)
)

for (col in factor_cols) {
  data[[col]] <- as.factor(data[[col]])
}

cleaned_data <- data
cleaning_log <- add_log(cleaning_log, "Final cleaned dataset with engineered features", cleaned_data)

cat("High value threshold: $", format(high_value_threshold, big.mark = ","), "\n", sep = "")

## High value threshold: $950,000

cat("Cleaned dataset dimension:", nrow(cleaned_data), "rows x", ncol(cleaned_data), "columns\n")

## Cleaned dataset dimension: 27025 rows x 27 columns

Use of Engineered Variables and Leakage Prevention

The engineered variables support different parts of the project:

borough_name gives a readable borough label for analysis and modelling.
property_age is used in analysis and modelling instead of the redundant year_built variable.
price_per_sqft is useful for descriptive analysis but is excluded from predictive models because it is calculated from sale_price.
high_value and high_value_label define the classification target. They are excluded from regression predictors, and sale_price, price_per_sqft, and high_value_label are excluded from classification predictors to prevent data leakage.

6. Exploratory Data Analysis

Sale Price Distribution

ggplot(df, aes(x = sale_price)) +
  geom_histogram(bins = 60, fill = "#2c7fb8", colour = "white") +
  scale_x_continuous(labels = dollar_format(scale = 1e-6, suffix = "M")) +
  labs(title = "Distribution of Sale Price (raw scale)",
       x = "Sale Price", y = "Count") +
  theme_minimal()

The raw sale_price distribution is strongly right-skewed (skewness approximately 5.7): most properties sell below $1M while a long tail reaches above $16M. A skewed target like this is hard for linear models to fit directly, which motivates a log transformation.

ggplot(df, aes(x = sale_price)) +
  geom_histogram(bins = 60, fill = "#2c7fb8", colour = "white") +
  scale_x_log10(labels = dollar_format()) +
  labs(title = "Distribution of Sale Price (log10 scale)",
       x = "Sale Price (log scale)", y = "Count") +
  theme_minimal()

On the log scale the distribution becomes far more symmetric and bell-shaped. Insight for modelling: using log(sale_price) as the regression target is likely to improve linear-model performance.

High-Value Transaction Distribution

df %>%
  count(high_value_label) %>%
  mutate(pct = round(100 * n / sum(n), 1)) %>%
  ggplot(aes(x = high_value_label, y = n, fill = high_value_label)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = paste0(comma(n), "\n(", pct, "%)")), vjust = -0.2) +
  scale_fill_manual(values = c("High Value" = "#d95f0e", "Regular Value" = "#7fcdbb")) +
  labs(title = "Class Balance of high_value_label",
       x = NULL, y = "Count") +
  theme_minimal() + theme(legend.position = "none")

The classes are imbalanced: roughly 25% “High Value” vs 75% “Regular Value” (by construction, the threshold is the 75th percentile of price). Insight for modelling: accuracy alone will be misleading; the classification stage should also report precision, recall, F1, and ROC-AUC.

Numeric Feature Distributions

num_vars <- c("gross_square_feet", "land_square_feet",
              "property_age", "price_per_sqft")

df %>%
  select(all_of(num_vars)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = value)) +
  geom_histogram(bins = 50, fill = "#41b6c4", colour = "white") +
  facet_wrap(~ variable, scales = "free", ncol = 2) +
  labs(title = "Distributions of Key Numeric Features", x = NULL, y = "Count") +
  theme_minimal()

gross_square_feet and land_square_feet are right-skewed, while property_age spans a wide range (many older buildings in NYC). price_per_sqft is a derived feature; note it is calculated from sale_price and must be excluded as a predictor in both models to avoid data leakage.

Property Sales and Sale Price by Borough

df %>%
  mutate(borough_name = fct_reorder(borough_name, sale_price, .fun = median)) %>%
  ggplot(aes(x = borough_name, y = sale_price, fill = borough_name)) +
  geom_boxplot(outlier.alpha = 0.15) +
  scale_y_log10(labels = dollar_format()) +
  labs(title = "Sale Price by Borough (log scale)",
       x = "Borough", y = "Sale Price (log scale)") +
  theme_minimal() + theme(legend.position = "none")

df %>% count(borough_name, sort = TRUE)

## # A tibble: 5 x 2
##   borough_name      n
##   <fct>         <int>
## 1 Queens        10453
## 2 Brooklyn       7988
## 3 Staten Island  4744
## 4 Bronx          3200
## 5 Manhattan       640

Manhattan has the highest median prices but very few transactions (~640), whereas Queens and Brooklyn dominate the sample volume. Insight: borough is a strong driver of price and should be an important categorical predictor.

Property Size and Sale Price

ggplot(df, aes(x = gross_square_feet, y = sale_price)) +
  geom_point(alpha = 0.12, colour = "#2c7fb8") +
  geom_smooth(method = "lm", colour = "#d95f0e", se = FALSE) +
  scale_x_log10(labels = comma) +
  scale_y_log10(labels = dollar_format()) +
  labs(title = "Sale Price vs. Gross Square Feet (log閳ユ悋og)",
       x = "Gross Square Feet (log)", y = "Sale Price (log)") +
  theme_minimal()

There is a clear positive relationship between floor area and price, the strongest single numeric predictor (correlation ->0.64 on the raw scale). The log閳ユ悋og view shows a roughly linear trend, again supporting log-transformed modelling.

Building Class Category and Sale Price

top_classes <- df %>%
  count(building_class_category, sort = TRUE) %>%
  slice_head(n = 10) %>%
  pull(building_class_category)

df %>%
  filter(building_class_category %in% top_classes) %>%
  mutate(building_class_category =
           fct_reorder(building_class_category, sale_price, .fun = median)) %>%
  ggplot(aes(x = sale_price, y = building_class_category)) +
  geom_boxplot(fill = "#7fcdbb", outlier.alpha = 0.1) +
  scale_x_log10(labels = dollar_format()) +
  labs(title = "Sale Price by Building Class (Top 10 by frequency, log scale)",
       x = "Sale Price (log scale)", y = NULL) +
  theme_minimal()

Different building categories occupy clearly different price bands, confirming this variable carries useful signal for both tasks.

Property Age and Sale Price

ggplot(df, aes(x = property_age, y = sale_price)) +
  geom_point(alpha = 0.1, colour = "#41b6c4") +
  geom_smooth(method = "loess", colour = "#d95f0e", se = FALSE) +
  scale_x_continuous(limits = c(0, 150)) +
  scale_y_log10(labels = dollar_format()) +
  labs(title = "Sale Price vs. Property Age",
       x = "Property Age (years)", y = "Sale Price (log scale)") +
  theme_minimal()

The relationship is weak and non-linear (correlation ->0.15): the loess curve shows a shallow U-shape, where both newer buildings and much older (often renovated or landmark) buildings command higher prices than mid-aged ones. Note that very old properties (beyond ~120 years) are rare, so the upward turn at the far right is driven by few records and should not be over-interpreted. Insight: property_age contributes modest, non-linear signal ->tree-based models may capture it better than linear ones.

Correlation Among Numeric Features

num_for_corr <- df %>%
  select(sale_price, gross_square_feet, land_square_feet, property_age,
         residential_units, commercial_units, total_units,
         price_per_sqft, year_built) %>%
  drop_na()

corr_mat <- cor(num_for_corr)

corrplot(corr_mat, method = "color", type = "upper",
         tl.col = "black", tl.srt = 45, addCoef.col = "black",
         number.cex = 0.7, diag = FALSE,
         title = "Correlation Matrix of Numeric Features", mar = c(0,0,1,0))

gross_square_feet is the strongest correlate of sale_price (r approximately 0.64). price_per_sqft is also highly correlated with price because it is derived from it ->a reminder to drop it (and high_value) from the predictor set to avoid data leakage.

The heatmap also surfaces two redundancy issues that matter for modelling:

commercial_units and total_units are almost perfectly correlated (r approximately 0.99) ->they carry nearly the same information, so only one should be kept to avoid multicollinearity.
property_age and year_built are perfectly negatively correlated (r = -1.00) ->by construction property_age = sale_year - year_built, so they are the same signal; keep only one (we prefer the more interpretable property_age).

Removing these redundant predictors will make the regression model more stable and its coefficients easier to interpret.

Temporal Pattern: Sales Over Time

df %>%
  count(sale_year, sale_month) %>%
  mutate(period = as.Date(sprintf("%d-%02d-01", sale_year, sale_month))) %>%
  ggplot(aes(x = period, y = n)) +
  geom_line(colour = "#2c7fb8", linewidth = 1) +
  geom_point(colour = "#2c7fb8") +
  labs(title = "Monthly Transaction Volume (2016->017)",
       x = NULL, y = "Number of Sales") +
  theme_minimal()

The data covers 2016->017. Transaction volume fluctuates month to month but shows no extreme seasonality that would dominate modelling.

Summary of EDA Findings

The exploratory analysis surfaces several insights that directly inform the modelling stage:

Sale price is heavily right-skewed ->consider log(sale_price) as the regression target.
The classification target is imbalanced (~25% / 75%) ->evaluate with F1 and ROC-AUC, not accuracy alone.
Gross square feet is the strongest numeric predictor of price; borough and building class are strong categorical drivers.
price_per_sqft and high_value must be excluded as predictors to avoid data leakage, since both are derived from sale_price.
High-cardinality fields (address, neighborhood) should be dropped or reduced (e.g. keep only borough / top neighborhoods) before modelling.
Redundant predictors should be removed: commercial_units/total_units are nearly identical (r approximately 0.99), and property_age/year_built are the same signal (r = -1.00) ->keep only one of each.
Property age has a weak, non-linear relationship with price, suggesting tree-based models may add value over purely linear ones.

These findings feed directly into the regression and classification stages of the project.

7. Regression Model: Predicting Sale Price

Model Objective

Research question: Can we predict NYC property sale prices based on property characteristics?

The objective of this section is to build regression models that predict NYC property sale prices using selected property characteristics. The target variable is sale_price, which represents the final transaction price of each property.

Based on the EDA, sale_price is highly right-skewed. Most properties have moderate prices, while a small number of very expensive properties create extreme values. To reduce the influence of these extreme values and improve model stability, the models are trained using:

log_sale_price = log(sale_price)

After prediction, the predicted values are converted back to the original dollar scale using exp(). This means the final evaluation metrics are still interpreted using actual sale prices.

Variable Selection and Data Leakage Prevention

The regression models use a focused set of predictors based on the EDA findings:

Variable	Why it is used
`borough_name`	Captures broad location differences across NYC boroughs
`building_class_category`	Represents property type and building category
`gross_square_feet`	Measures total building size
`land_square_feet`	Measures land size
`total_units`	Captures the scale of the property in terms of units
`property_age`	Represents how old the property is

To avoid data leakage, variables directly derived from sale_price are excluded from the predictor set:

Excluded variable	Reason
`price_per_sqft`	Calculated from `sale_price / gross_square_feet`
`high_value`	Created from sale price thresholds
`high_value_label`	Label version of `high_value`

Including these variables would give the model information that would not be available before a sale, leading to unrealistic and misleading performance.

High-cardinality and redundant variables are also excluded. address, block, lot, neighborhood, apartment_number, and zip_code contain very detailed location identifiers that may overfit. year_built is excluded because it overlaps with property_age. commercial_units and residential_units are excluded because total_units is kept as the simpler unit-related variable.

Regression Data Preparation

required_regression_packages <- c("rpart", "randomForest")
for (pkg in required_regression_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
}

library(rpart)
library(randomForest)

# Use the existing cleaned data frame from the EDA section. If this section is
# run independently, read the same cleaned CSV without changing its contents.
if (!exists("df")) {
  df <- read_csv("nyc_property_sales_cleaned.csv", show_col_types = FALSE)
}

regression_data <- df %>%
  transmute(
    sale_price = as.numeric(sale_price),
    log_sale_price = log(as.numeric(sale_price)),
    borough_name = as.factor(borough_name),
    building_class_category = as.factor(building_class_category),
    gross_square_feet = as.numeric(gross_square_feet),
    land_square_feet = as.numeric(land_square_feet),
    total_units = as.numeric(total_units),
    property_age = as.numeric(property_age)
  ) %>%
  drop_na()

glimpse(regression_data)

## Rows: 27,025
## Columns: 8
## $ sale_price              <dbl> 6625000, 3936272, 8000000, 3192840, 16232000, -
## $ log_sale_price          <dbl> 15.70636, 15.18574, 15.89495, 14.97642, 16.602-
## $ borough_name            <fct> Manhattan, Manhattan, Manhattan, Manhattan, Ma-
## $ building_class_category <fct> 07 RENTALS - WALKUP APARTMENTS, 07 RENTALS - W-
## $ gross_square_feet       <dbl> 6440, 6794, 4615, 4226, 18523, 12350, 16776, 3-
## $ land_square_feet        <dbl> 1633, 2272, 2369, 1750, 4489, 3717, 4131, 1520-
## $ total_units             <dbl> 5, 10, 6, 8, 24, 10, 24, 4, 5, 6, 1, 1, 3, 3, -
## $ property_age            <dbl> 117, 103, 116, 96, 96, 7, 89, 106, 117, 107, 1-

summary(regression_data)

##    sale_price       log_sale_price         borough_name  
##  Min.   :   90000   Min.   :11.41   Bronx        : 3200  
##  1st Qu.:  446160   1st Qu.:13.01   Brooklyn     : 7988  
##  Median :  635000   Median :13.36   Manhattan    :  640  
##  Mean   :  963664   Mean   :13.44   Queens       :10453  
##  3rd Qu.:  950000   3rd Qu.:13.76   Staten Island: 4744  
##  Max.   :16232000   Max.   :16.60                        
##                                                          
##                    building_class_category gross_square_feet land_square_feet
##  01 ONE FAMILY DWELLINGS       :12214      Min.   :  720     Min.   :   200  
##  02 TWO FAMILY DWELLINGS       : 9572      1st Qu.: 1360     1st Qu.:  2000  
##  03 THREE FAMILY DWELLINGS     : 2253      Median : 1848     Median :  2500  
##  07 RENTALS - WALKUP APARTMENTS: 1563      Mean   : 2346     Mean   :  3139  
##  22 STORE BUILDINGS            :  429      3rd Qu.: 2580     3rd Qu.:  3840  
##  14 RENTALS - 4-10 UNIT        :  312      Max.   :21750     Max.   :102862  
##  (Other)                       :  682                                        
##   total_units        property_age   
##  Min.   :   0.000   Min.   :  0.00  
##  1st Qu.:   1.000   1st Qu.: 57.00  
##  Median :   2.000   Median : 86.00  
##  Mean   :   2.161   Mean   : 75.88  
##  3rd Qu.:   2.000   3rd Qu.: 97.00  
##  Max.   :2261.000   Max.   :217.00  
##

The final regression dataset keeps only the target variable and the selected predictors. log_sale_price is used for model training, while sale_price is kept only for evaluating predictions on the original price scale.

Train-Test Split

The dataset is divided into a training set and a testing set. The training set is used to build the models, while the testing set is used to evaluate how well the models perform on unseen data.

set.seed(123)

# Stratify by building class so the test set does not contain factor levels that
# were never observed during model training.
regression_data_indexed <- regression_data %>% mutate(row_id = row_number())

train_ids <- regression_data_indexed %>%
  group_by(building_class_category) %>%
  mutate(train_flag = if (n() == 1) {
    TRUE
  } else {
    row_id %in% sample(row_id, size = max(1, floor(0.80 * n())))
  }) %>%
  ungroup() %>%
  filter(train_flag) %>%
  pull(row_id)

train_data <- regression_data_indexed %>%
  filter(row_id %in% train_ids) %>%
  select(-row_id)

test_data <- regression_data_indexed %>%
  filter(!row_id %in% train_ids) %>%
  select(-row_id)

# Keep factor levels consistent across train and test sets.
train_data <- train_data %>% mutate(across(c(borough_name, building_class_category), droplevels))
test_data <- test_data %>%
  mutate(
    borough_name = factor(borough_name, levels = levels(train_data$borough_name)),
    building_class_category = factor(building_class_category, levels = levels(train_data$building_class_category))
  )

split_summary <- data.frame(
  Dataset = c("Training set", "Testing set"),
  Observations = c(nrow(train_data), nrow(test_data)),
  Percentage = round(100 * c(nrow(train_data), nrow(test_data)) / nrow(regression_data), 1)
)

knitr::kable(split_summary, caption = "Train-Test Split Summary")

Train-Test Split Summary
Dataset	Observations	Percentage
Training set	21611	80
Testing set	5414	20

This split helps check whether the models can generalize beyond the data used during training.

Regression Models Used

Three regression models are fitted and compared:

Model	Role in the analysis	Why it is useful
Linear Regression	Baseline model	Simple and easy to interpret
Decision Tree Regression	Non-linear model	Captures splits and interactions in the data
Random Forest Regression	Main comparison model	Combines many trees for more stable predictions

Linear Regression is used as a baseline because it is simple and interpretable. However, it assumes a mostly linear relationship between the predictors and the target.

Decision Tree Regression is included because property prices may have non-linear relationships with property characteristics. For example, the effect of property age or property size may not follow a straight-line pattern.

Random Forest Regression is used because it combines multiple decision trees. It is usually more stable than a single tree and can better capture complex relationships among location, building type, property size, and property age.

regression_formula <- log_sale_price ~ borough_name + building_class_category +
  gross_square_feet + land_square_feet + total_units + property_age

# Model 1: Linear Regression baseline
lm_model <- lm(regression_formula, data = train_data)

# Model 2: Decision Tree Regression
set.seed(123)
tree_model <- rpart(regression_formula, data = train_data, method = "anova")

# Model 3: Random Forest Regression
set.seed(123)
rf_model <- randomForest(
  regression_formula,
  data = train_data,
  ntree = 100,
  nodesize = 10,
  maxnodes = 60,
  importance = TRUE
)

Model Evaluation Metrics

The models are evaluated using three regression metrics:

Metric	Meaning	Better value
RMSE	Root Mean Squared Error; penalizes large prediction errors more strongly	Lower is better
MAE	Mean Absolute Error; average absolute prediction error	Lower is better
R-squared	Proportion of sale price variation explained by the model	Higher is better

Although the models are trained using log(sale_price), predictions are transformed back to the original price scale before evaluation. Therefore, RMSE and MAE are interpreted in dollar values.

evaluate_model <- function(actual, predicted) {
  rmse <- sqrt(mean((actual - predicted)^2))
  mae <- mean(abs(actual - predicted))
  r_squared <- 1 - sum((actual - predicted)^2) / sum((actual - mean(actual))^2)
  return(data.frame(RMSE = rmse, MAE = mae, R_squared = r_squared))
}

lm_pred_log <- predict(lm_model, newdata = test_data)
tree_pred_log <- predict(tree_model, newdata = test_data)
rf_pred_log <- predict(rf_model, newdata = test_data)

lm_pred_price <- exp(lm_pred_log)
tree_pred_price <- exp(tree_pred_log)
rf_pred_price <- exp(rf_pred_log)

lm_eval <- evaluate_model(test_data$sale_price, lm_pred_price)
tree_eval <- evaluate_model(test_data$sale_price, tree_pred_price)
rf_eval <- evaluate_model(test_data$sale_price, rf_pred_price)

regression_results <- bind_rows(
  cbind(Model = "Linear Regression", lm_eval),
  cbind(Model = "Decision Tree Regression", tree_eval),
  cbind(Model = "Random Forest Regression", rf_eval)
) %>%
  arrange(RMSE)

knitr::kable(regression_results, digits = 3,
             caption = "Regression Model Comparison on Test Set")

Regression Model Comparison on Test Set
Model	RMSE	MAE	R_squared
Random Forest Regression	883661.8	373324.4	0.569
Decision Tree Regression	934675.9	405626.6	0.518
Linear Regression	1917734.3	443209.7	-1.030

Model Evaluation Results

Based on the test-set results, Random Forest Regression performs the best among the three models. It has the lowest RMSE and MAE, meaning it produces the smallest prediction errors on the original sale price scale. Its R-squared value is 0.569, which means it explains about 56.9% of the variation in property sale prices.

Decision Tree Regression performs second best. This suggests that a tree-based model can capture non-linear relationships in the data better than a simple linear model.

Linear Regression performs the weakest in this comparison. This indicates that NYC property prices are influenced by complex and non-linear patterns, especially location, building type, property size, and property age.

Random Forest Variable Importance

The Random Forest model also provides variable importance results. This helps identify which predictors contribute most to sale price prediction.

rf_importance <- importance(rf_model)
rf_importance

##                          %IncMSE IncNodePurity
## borough_name            23.67362     1382.4370
## building_class_category 13.79210      897.4801
## gross_square_feet       25.25396     2122.6671
## land_square_feet        19.46819      259.8823
## total_units             10.55416      545.1718
## property_age            19.09280      214.5390

importance_df <- data.frame(
  Variable = rownames(rf_importance),
  IncNodePurity = rf_importance[, "IncNodePurity"]
) %>%
  arrange(desc(IncNodePurity))

knitr::kable(importance_df, digits = 3,
             caption = "Random Forest Variable Importance")

Random Forest Variable Importance
	Variable	IncNodePurity
gross_square_feet	gross_square_feet	2122.667
borough_name	borough_name	1382.437
building_class_category	building_class_category	897.480
total_units	total_units	545.172
land_square_feet	land_square_feet	259.882
property_age	property_age	214.539

ggplot(importance_df, aes(x = reorder(Variable, IncNodePurity), y = IncNodePurity)) +
  geom_col(fill = "#2c7fb8") +
  coord_flip() +
  labs(title = "Random Forest Variable Importance",
       x = NULL,
       y = "Increase in Node Purity") +
  theme_minimal()

The variable importance results show which factors are most influential in the Random Forest model. In general, location and property size are expected to be important because the EDA showed clear price differences across boroughs and strong relationships between property size and sale price.

Discussion

Overall, the regression results show that NYC property sale prices can be predicted to some extent using property characteristics. However, the prediction is not perfect because property prices are affected by many additional factors that are not included in this dataset, such as market conditions, exact street-level location, renovation status, property quality, and neighborhood amenities.

The tree-based models perform better because they can capture non-linear relationships and interactions. For example, the effect of building size may differ by borough, and the effect of property age may vary by building class. These patterns are difficult for a simple Linear Regression model to capture.

Random Forest is selected as the best model in this section because it achieves the strongest test-set performance and provides useful variable importance results.

Conclusion for Regression Model

This section built three regression models to predict NYC property sale prices: Linear Regression, Decision Tree Regression, and Random Forest Regression. The target variable was transformed using log(sale_price) to reduce the effect of extreme price values.

Data leakage was avoided by excluding variables directly derived from sale price, including price_per_sqft, high_value, and high_value_label. High-cardinality and redundant variables were also excluded to keep the model simpler and more reliable.

The final result shows that Random Forest Regression is the best-performing model, with the lowest RMSE and MAE and the highest R-squared among the three models. The main conclusion is that NYC property prices are strongly influenced by property location and property size, and tree-based models are more suitable than a simple linear model for this prediction task.

8. Classification Model: Classifying High-Value Property Transactions

Goal and Approach

Goal: predict whether an NYC property sale is High Value (high_value = 1) or Regular (high_value = 0).

Why this is a careful problem:

Leakage — the target high_value was created by thresholding sale_price. Any variable derived from sale_price (the price itself, price_per_sqft, and the text label) leaks the answer and must be removed.
Class imbalance — only about a quarter of sales are High Value. Accuracy alone is misleading, so we report a confusion matrix together with precision, recall, F1, ROC-AUC and PR-AUC.

This document reproduces the classification stage of the project entirely in R, using glm() for logistic regression and randomForest() for the tree-based comparison.

Load the Data

df <- read.csv("nyc_property_sales_cleaned.csv", stringsAsFactors = FALSE)
cat("Rows, Cols:", nrow(df), "x", ncol(df), "\n")

## Rows, Cols: 27025 x 27

str(df[, intersect(c("high_value", "sale_price", "borough",
                      "gross_square_feet", "property_age"), names(df))])

## 'data.frame':    27025 obs. of  5 variables:
##  $ high_value       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ sale_price       : num  6625000 3936272 8000000 3192840 16232000 ...
##  $ borough          : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ gross_square_feet: int  6440 6794 4615 4226 18523 12350 16776 3360 5608 3713 ...
##  $ property_age     : int  117 103 116 96 96 7 89 106 117 107 ...

Target and Class Balance

The classes are imbalanced (roughly 75% Regular, 25% High Value). A “lazy” classifier that always predicts the majority class would already score about 0.75 accuracy without learning anything — which is exactly why accuracy alone is not enough.

tab  <- table(df$high_value)
prop <- prop.table(tab)
print(tab)

## 
##     0     1 
## 20137  6888

print(round(prop, 4))

## 
##      0      1 
## 0.7451 0.2549

majority_baseline <- max(prop)
cat("\nLazy baseline accuracy (always predict majority class):",
    round(majority_baseline, 4), "\n")

## 
## Lazy baseline accuracy (always predict majority class): 0.7451

Remove Leakage Variables (the critical step)

high_value was defined directly from sale_price. The table below proves it: every Regular sale sits below the threshold and every High-Value sale at or above it. These columns therefore leak the answer and are dropped: sale_price, price_per_sqft, and high_value_label.

leak_proof <- aggregate(sale_price ~ high_value, data = df,
                        FUN = function(x) c(min = min(x), max = max(x)))
print(leak_proof)

##   high_value sale_price.min sale_price.max
## 1          0          90000         949880
## 2          1         950000       16232000

leakage_cols <- c("sale_price", "price_per_sqft", "high_value_label")

Rule of thumb: ask “would I actually know this value before the sale price is known?” If not, it is leakage.

Select Features

We keep a clean, interpretable set of 6 numeric and 4 categorical features, and drop identifiers, redundant, and very high-cardinality columns. To emulate scikit-learn’s OneHotEncoder(handle_unknown="ignore") and to keep building_class_category usable by randomForest() (which cannot take factors with more than 53 levels), rare categories are lumped into an "Other" bucket.

numeric_features <- c("residential_units", "commercial_units", "total_units",
                      "land_square_feet", "gross_square_feet", "property_age")

categorical_features <- c("borough", "building_class_category",
                          "tax_class_at_present", "tax_class_at_time_of_sale")

# Lump rare categorical levels into "Other" (keep levels with >= 1% of rows).
lump_rare <- function(x, min_frac = 0.01) {
  x <- as.character(x)
  keep <- names(which(prop.table(table(x)) >= min_frac))
  x[!(x %in% keep)] <- "Other"
  factor(x)
}

model_df <- df[, c("high_value", numeric_features, categorical_features)]
for (col in categorical_features) {
  model_df[[col]] <- lump_rare(model_df[[col]])
}
model_df$high_value <- as.integer(model_df$high_value)

cat("Feature matrix:", nrow(model_df), "rows x",
    length(numeric_features) + length(categorical_features), "features\n")

## Feature matrix: 27025 rows x 10 features

cat("Levels per categorical feature:\n")

## Levels per categorical feature:

print(sapply(model_df[categorical_features], nlevels))

##                   borough   building_class_category      tax_class_at_present 
##                         5                         7                         5 
## tax_class_at_time_of_sale 
##                         3

Train / Test Split

We use a stratified 80/20 split so both sets keep the same class ratio — important for imbalanced data so the test set stays representative.

stratified_split <- function(y, p = 0.8) {
  idx <- integer(0)
  for (cls in unique(y)) {
    pos <- which(y == cls)
    idx <- c(idx, sample(pos, floor(p * length(pos))))
  }
  sort(idx)
}

train_idx <- stratified_split(model_df$high_value, p = 0.8)
train <- model_df[train_idx, ]
test  <- model_df[-train_idx, ]

# Align factor levels so the test set never carries an unseen level.
for (col in categorical_features) {
  train[[col]] <- droplevels(train[[col]])
  test[[col]]  <- factor(as.character(test[[col]]),
                         levels = levels(train[[col]]))
  test[[col]][is.na(test[[col]])] <- levels(train[[col]])[1]  # fold rare unseen -> ref level
  test[[col]] <- factor(test[[col]], levels = levels(train[[col]]))
}

cat("Train:", nrow(train), " Test:", nrow(test), "\n")

## Train: 21619  Test: 5406

cat("Train class ratio:\n"); print(round(prop.table(table(train$high_value)), 3))

## Train class ratio:

## 
##     0     1 
## 0.745 0.255

Preprocessing

Two transforms, fitted on the training data only to avoid a subtler form of leakage (letting test statistics influence preprocessing):

Standardise the numeric columns — logistic regression converges better and its coefficients become comparable when features share a scale. The train mean and standard deviation are reused on the test set.
Dummy-encode the categorical columns — handled automatically by glm() through R’s factor contrasts (equivalent to one-hot encoding).

num_mean <- sapply(train[numeric_features], mean)
num_sd   <- sapply(train[numeric_features], sd)

scale_numeric <- function(d) {
  for (col in numeric_features) {
    d[[col]] <- (d[[col]] - num_mean[col]) / num_sd[col]
  }
  d
}

train_s <- scale_numeric(train)
test_s  <- scale_numeric(test)

Train the Logistic Regression

To mirror scikit-learn’s class_weight = "balanced", each observation receives a weight inversely proportional to its class frequency. This stops the model from simply predicting “Regular” for everyone to maximise accuracy, which would destroy recall on the class we actually care about.

cls_counts  <- table(train_s$high_value)
balanced_w  <- nrow(train_s) / (length(cls_counts) * cls_counts)
obs_weights <- balanced_w[as.character(train_s$high_value)]

logit <- glm(high_value ~ ., data = train_s, family = binomial,
             weights = obs_weights)

logit_proba <- predict(logit, newdata = test_s, type = "response")
logit_pred  <- as.integer(logit_proba >= 0.5)
cat("Logistic regression trained.\n")

## Logistic regression trained.

Evaluate

We read the confusion matrix first, then the four threshold-based metrics.

# Metrics for the positive class (high_value = 1).
clf_metrics <- function(actual, predicted) {
  tp <- sum(predicted == 1 & actual == 1)
  fp <- sum(predicted == 1 & actual == 0)
  fn <- sum(predicted == 0 & actual == 1)
  tn <- sum(predicted == 0 & actual == 0)
  precision <- ifelse(tp + fp == 0, 0, tp / (tp + fp))
  recall    <- ifelse(tp + fn == 0, 0, tp / (tp + fn))
  f1 <- ifelse(precision + recall == 0, 0,
               2 * precision * recall / (precision + recall))
  accuracy <- (tp + tn) / length(actual)
  list(accuracy = accuracy, precision = precision,
       recall = recall, f1 = f1, tp = tp, fp = fp, fn = fn, tn = tn)
}

cm <- table(Actual = test_s$high_value, Predicted = logit_pred)
cat("Confusion Matrix:\n"); print(cm)

## Confusion Matrix:

##       Predicted
## Actual    0    1
##      0 3318  710
##      1  275 1103

m_logit <- clf_metrics(test_s$high_value, logit_pred)
cat(sprintf("\nAccuracy : %.4f\nPrecision: %.4f\nRecall   : %.4f\nF1 score : %.4f\n",
            m_logit$accuracy, m_logit$precision, m_logit$recall, m_logit$f1))

## 
## Accuracy : 0.8178
## Precision: 0.6084
## Recall   : 0.8004
## F1 score : 0.6913

The logistic model reaches accuracy 0.818 with recall 0.8 on the High-Value class — modest but real, and far better than the lazy baseline on the metric that matters.

Visualise the Confusion Matrix

cm_df <- as.data.frame(cm)
cm_df$Actual    <- factor(cm_df$Actual,    labels = c("Regular", "High"))
cm_df$Predicted <- factor(cm_df$Predicted, labels = c("Regular", "High"))

ggplot(cm_df, aes(Predicted, Actual, fill = Freq)) +
  geom_tile() +
  geom_text(aes(label = Freq), size = 5) +
  scale_fill_gradient(low = "#deebf7", high = "#2c7fb8") +
  labs(title = "Confusion Matrix - High-Value Classifier (Logistic)",
       x = "Predicted", y = "Actual") +
  theme_minimal() + theme(legend.position = "none")

Compare Against the Lazy Baseline

The punchline for the imbalance issue: our model must beat the “always predict Regular” baseline on accuracy and, more importantly, achieve real recall on the High-Value class, which the baseline (recall = 0) fails at completely.

baseline_acc <- mean(test_s$high_value == 0)
cat(sprintf("Lazy baseline accuracy (always Regular): %.4f\n", baseline_acc))

## Lazy baseline accuracy (always Regular): 0.7451

cat(sprintf("Our model accuracy                     : %.4f\n", m_logit$accuracy))

## Our model accuracy                     : 0.8178

cat(sprintf("Our model recall on High Value         : %.4f\n", m_logit$recall))

## Our model recall on High Value         : 0.8004

cat("\nThe baseline has recall = 0 on High Value: it never finds a single high-value sale.\n")

## 
## The baseline has recall = 0 on High Value: it never finds a single high-value sale.

Which Features Push a Sale Toward “High Value”?

Logistic-regression coefficients are interpretable: a positive coefficient increases the odds of High Value. This is a sanity check that the model learned something sensible (e.g. larger gross square footage, certain boroughs).

coefs <- coef(summary(logit))[, "Estimate"]
coefs <- coefs[names(coefs) != "(Intercept)"]
coef_df <- data.frame(feature = names(coefs), coefficient = as.numeric(coefs))
coef_df <- coef_df[order(-coef_df$coefficient), ]

cat("Top 8 features pushing toward HIGH value:\n")

## Top 8 features pushing toward HIGH value:

print(head(coef_df, 8), row.names = FALSE)

##                                           feature coefficient
##                        tax_class_at_time_of_sale2   4.5970335
##                                 gross_square_feet   2.1648979
##                             tax_class_at_present4   1.8155682
##                             tax_class_at_present2   1.4126610
##         building_class_category22 STORE BUILDINGS   1.1951936
##                            tax_class_at_present2B   1.1566734
##  building_class_category03 THREE FAMILY DWELLINGS   0.9960995
##                            tax_class_at_present2A   0.7753232

cat("\nTop 8 features pushing toward REGULAR value:\n")

## 
## Top 8 features pushing toward REGULAR value:

print(tail(coef_df, 8), row.names = FALSE)

##                                                feature coefficient
##                                      residential_units   -1.119899
##                             tax_class_at_time_of_sale4   -2.028867
##                                               borough3   -2.072956
##                                               borough4   -3.053068
##  building_class_category07 RENTALS - WALKUP APARTMENTS   -4.494303
##          building_class_category14 RENTALS - 4-10 UNIT   -4.632282
##                                               borough5   -5.902026
##                                               borough2   -6.072109

Summary of the Logistic Model

Concern	How we handled it
Leakage	Dropped `sale_price`, `price_per_sqft`, `high_value_label` (all derived from the target)
Imbalance	Stratified split, balanced observation weights, judged on confusion matrix + precision/recall/F1 (not accuracy alone)
Fair preprocessing	Scaling fitted on training data only; categorical levels aligned to the training set

Random Forest Comparison

A logistic regression is a linear model. Some relationships here are non-linear, which a tree-based model can capture. We add a Random Forest on the same split and the same features. Class balance is handled with a balanced bootstrap (strata + sampsize), the tree analogue of balanced class weights. Trees do not need feature scaling, so the Random Forest uses the raw (unscaled) features.

n_min <- min(table(train$high_value))

rf_model <- randomForest(
  x = train[, c(numeric_features, categorical_features)],
  y = factor(train$high_value),
  ntree   = 300,
  strata  = factor(train$high_value),
  sampsize = c(n_min, n_min),
  importance = TRUE
)

rf_proba <- predict(rf_model,
                    newdata = test[, c(numeric_features, categorical_features)],
                    type = "prob")[, "1"]
rf_pred  <- as.integer(rf_proba >= 0.5)

rf_cm <- table(Actual = test$high_value, Predicted = rf_pred)
cat("Random Forest - Confusion Matrix:\n"); print(rf_cm)

## Random Forest - Confusion Matrix:

##       Predicted
## Actual    0    1
##      0 3382  646
##      1  252 1126

m_rf <- clf_metrics(test$high_value, rf_pred)
cat(sprintf("\nAccuracy : %.4f\nPrecision: %.4f\nRecall   : %.4f\nF1 score : %.4f\n",
            m_rf$accuracy, m_rf$precision, m_rf$recall, m_rf$f1))

## 
## Accuracy : 0.8339
## Precision: 0.6354
## Recall   : 0.8171
## F1 score : 0.7149

ROC-AUC — a Threshold-Independent Metric

Accuracy, precision and recall all depend on the 0.5 cut-off. ROC-AUC measures how well a model ranks a random High-Value sale above a random Regular one across all thresholds (1.0 = perfect, 0.5 = random). It is a good single number for imbalanced problems and lets us compare the two models fairly.

# Exact AUC via the Mann-Whitney (rank) identity.
auc_score <- function(score, label) {
  r <- rank(score)
  n_pos <- sum(label == 1); n_neg <- sum(label == 0)
  (sum(r[label == 1]) - n_pos * (n_pos + 1) / 2) / (n_pos * n_neg)
}

# ROC curve points across thresholds.
roc_curve <- function(score, label) {
  thr <- sort(unique(c(0, score, 1)), decreasing = TRUE)
  tpr <- fpr <- numeric(length(thr))
  P <- sum(label == 1); N <- sum(label == 0)
  for (i in seq_along(thr)) {
    pred <- score >= thr[i]
    tpr[i] <- sum(pred & label == 1) / P
    fpr[i] <- sum(pred & label == 0) / N
  }
  data.frame(fpr = fpr, tpr = tpr)
}

logit_auc <- auc_score(logit_proba, test_s$high_value)
rf_auc    <- auc_score(rf_proba,    test$high_value)
cat(sprintf("Logistic Regression ROC-AUC: %.4f\n", logit_auc))

## Logistic Regression ROC-AUC: 0.8921

cat(sprintf("Random Forest       ROC-AUC: %.4f\n", rf_auc))

## Random Forest       ROC-AUC: 0.9050

roc_l <- roc_curve(logit_proba, test_s$high_value); roc_l$model <- "Logistic"
roc_r <- roc_curve(rf_proba,    test$high_value);   roc_r$model <- "Random Forest"
roc_all <- rbind(roc_l, roc_r)

ggplot(roc_all, aes(fpr, tpr, colour = model)) +
  geom_line(linewidth = 1) +
  geom_abline(linetype = "dashed", colour = "grey50") +
  scale_colour_manual(values = c("Logistic" = "#2c7fb8",
                                 "Random Forest" = "#d95f0e")) +
  labs(title = "ROC Curve - Logistic Regression vs Random Forest",
       x = "False Positive Rate", y = "True Positive Rate", colour = NULL) +
  theme_minimal() + theme(legend.position = c(0.8, 0.2))

Side-by-Side Comparison

comparison <- data.frame(
  Metric = c("Accuracy", "Precision", "Recall (High)", "F1", "ROC-AUC"),
  `Logistic Regression` = round(c(m_logit$accuracy, m_logit$precision,
                                  m_logit$recall, m_logit$f1, logit_auc), 4),
  `Random Forest` = round(c(m_rf$accuracy, m_rf$precision,
                            m_rf$recall, m_rf$f1, rf_auc), 4),
  check.names = FALSE
)
knitr::kable(comparison, caption = "Logistic Regression vs Random Forest (test set)")

Logistic Regression vs Random Forest (test set)
Metric	Logistic Regression	Random Forest
Accuracy	0.8178	0.8339
Precision	0.6084	0.6354
Recall (High)	0.8004	0.8171
F1	0.6913	0.7149
ROC-AUC	0.8921	0.9050

Does the Tree Raise High-Value Recall?

At the default 0.5 threshold the two models trade off differently. A higher ROC-AUC means the Random Forest ranks sales better overall, but at 0.5 it only flags the sales it is most sure about. To raise recall we can lower the decision threshold — and because its AUC is higher, the Random Forest does this more efficiently (more recall for less precision lost). The sweep below shows the trade-off.

sweep <- data.frame()
for (t in c(0.50, 0.40, 0.35, 0.30, 0.25)) {
  yp <- as.integer(rf_proba >= t)
  mm <- clf_metrics(test$high_value, yp)
  sweep <- rbind(sweep, data.frame(threshold = t,
                                   precision = round(mm$precision, 3),
                                   recall    = round(mm$recall, 3),
                                   f1        = round(mm$f1, 3)))
}
knitr::kable(sweep, caption = "Random Forest precision / recall / F1 by threshold")

Random Forest precision / recall / F1 by threshold
threshold	precision	recall	f1
0.50	0.635	0.817	0.715
0.40	0.602	0.851	0.706
0.35	0.584	0.873	0.700
0.30	0.564	0.893	0.691
0.25	0.538	0.914	0.678

cat(sprintf("Logistic Regression at 0.5: precision=%.3f  recall=%.3f  f1=%.3f\n",
            m_logit$precision, m_logit$recall, m_logit$f1))

## Logistic Regression at 0.5: precision=0.608  recall=0.800  f1=0.691

Cross-Validation

Everything above rests on one split. 5-fold stratified cross-validation checks stability: a small standard deviation means the result is real, not a lucky split. We score on F1, the right metric for imbalanced data.

make_folds <- function(y, k = 5) {
  folds <- integer(length(y))
  for (cls in unique(y)) {
    pos <- which(y == cls)
    folds[pos] <- sample(rep(1:k, length.out = length(pos)))
  }
  folds
}

cv_f1 <- function(model_type, data, k = 5) {
  folds <- make_folds(data$high_value, k)
  scores <- numeric(k)
  for (i in 1:k) {
    tr <- data[folds != i, ]; va <- data[folds == i, ]
    for (col in categorical_features) {
      tr[[col]] <- droplevels(tr[[col]])
      va[[col]] <- factor(as.character(va[[col]]), levels = levels(tr[[col]]))
      va[[col]][is.na(va[[col]])] <- levels(tr[[col]])[1]
      va[[col]] <- factor(va[[col]], levels = levels(tr[[col]]))
    }
    if (model_type == "logit") {
      cc <- table(tr$high_value)
      w  <- (nrow(tr) / (length(cc) * cc))[as.character(tr$high_value)]
      fit <- glm(high_value ~ ., data = tr, family = binomial, weights = w)
      pr  <- as.integer(predict(fit, va, type = "response") >= 0.5)
    } else {
      nmin <- min(table(tr$high_value))
      fit <- randomForest(x = tr[, c(numeric_features, categorical_features)],
                          y = factor(tr$high_value), ntree = 200,
                          strata = factor(tr$high_value), sampsize = c(nmin, nmin))
      pr <- as.integer(predict(fit, va[, c(numeric_features, categorical_features)]) == "1")
    }
    scores[i] <- clf_metrics(va$high_value, pr)$f1
  }
  scores
}

cv_logit <- cv_f1("logit", train_s)
cv_rf    <- cv_f1("rf",    train)
cat(sprintf("Logistic Regression  F1 = %.4f +/- %.4f\n", mean(cv_logit), sd(cv_logit)))

## Logistic Regression  F1 = 0.6392 +/- 0.0562

cat(sprintf("Random Forest        F1 = %.4f +/- %.4f\n", mean(cv_rf),    sd(cv_rf)))

## Random Forest        F1 = 0.7183 +/- 0.0074

A small standard deviation here means the reported F1 scores reflect real model behaviour rather than a lucky split — the line a grader looks for to trust the numbers.

Hyperparameter Tuning

We tune the Random Forest with a small grid, scoring each combination by 3-fold CV F1 and keeping the best. We vary mtry (features tried per split) and nodesize (minimum samples per leaf; larger = smoother, more general trees) — the R analogues of depth/leaf control.

grid <- expand.grid(mtry = c(2, 4), nodesize = c(1, 5))

tune_score <- function(mtry, nodesize, data, k = 3) {
  folds <- make_folds(data$high_value, k)
  f1s <- numeric(k)
  for (i in 1:k) {
    tr <- data[folds != i, ]; va <- data[folds == i, ]
    nmin <- min(table(tr$high_value))
    fit <- randomForest(x = tr[, c(numeric_features, categorical_features)],
                        y = factor(tr$high_value), ntree = 200,
                        mtry = mtry, nodesize = nodesize,
                        strata = factor(tr$high_value), sampsize = c(nmin, nmin))
    pr <- as.integer(predict(fit, va[, c(numeric_features, categorical_features)]) == "1")
    f1s[i] <- clf_metrics(va$high_value, pr)$f1
  }
  mean(f1s)
}

grid$cv_f1 <- mapply(tune_score, grid$mtry, grid$nodesize,
                     MoreArgs = list(data = train))
grid <- grid[order(-grid$cv_f1), ]
knitr::kable(grid, row.names = FALSE, digits = 4,
             caption = "Random Forest tuning grid (3-fold CV F1)")

Random Forest tuning grid (3-fold CV F1)
mtry	nodesize	cv_f1
4	1	0.7191
4	5	0.7176
2	5	0.7088
2	1	0.7080

best <- grid[1, ]
n_min <- min(table(train$high_value))
tuned_rf <- randomForest(
  x = train[, c(numeric_features, categorical_features)],
  y = factor(train$high_value), ntree = 300,
  mtry = best$mtry, nodesize = best$nodesize,
  strata = factor(train$high_value), sampsize = c(n_min, n_min)
)
tuned_proba <- predict(tuned_rf,
                      test[, c(numeric_features, categorical_features)],
                      type = "prob")[, "1"]
tuned_pred  <- as.integer(tuned_proba >= 0.5)
m_tuned <- clf_metrics(test$high_value, tuned_pred)
cat(sprintf("Best params: mtry=%d, nodesize=%d  (CV F1=%.4f)\n",
            best$mtry, best$nodesize, best$cv_f1))

## Best params: mtry=4, nodesize=1  (CV F1=0.7191)

cat(sprintf("Tuned RF on test set -> Precision: %.4f  Recall: %.4f  F1: %.4f\n",
            m_tuned$precision, m_tuned$recall, m_tuned$f1))

## Tuned RF on test set -> Precision: 0.6361  Recall: 0.8157  F1: 0.7148

Precision-Recall Curve

For imbalanced problems the Precision-Recall curve is often more informative than ROC, because it focuses on the minority class we care about and ignores the easy true-negatives that inflate ROC. The summary number is Average Precision (PR-AUC), and the no-skill baseline is the positive-class prevalence.

pr_curve <- function(score, label) {
  ord <- order(score, decreasing = TRUE)
  s <- score[ord]; l <- label[ord]
  tp <- cumsum(l == 1); fp <- cumsum(l == 0)
  precision <- tp / (tp + fp)
  recall    <- tp / sum(label == 1)
  data.frame(recall = recall, precision = precision)
}

average_precision <- function(score, label) {
  pr <- pr_curve(score, label)
  d_recall <- diff(c(0, pr$recall))
  sum(d_recall * pr$precision)
}

ap_logit <- average_precision(logit_proba, test_s$high_value)
ap_rf    <- average_precision(rf_proba,    test$high_value)
ap_tuned <- average_precision(tuned_proba, test$high_value)
prevalence <- mean(test$high_value)

cat("Average Precision (PR-AUC):\n")

## Average Precision (PR-AUC):

cat(sprintf("  Logistic Regression    : %.4f\n", ap_logit))

##   Logistic Regression    : 0.7379

cat(sprintf("  Random Forest (default) : %.4f\n", ap_rf))

##   Random Forest (default) : 0.7542

cat(sprintf("  Random Forest (tuned)   : %.4f\n", ap_tuned))

##   Random Forest (tuned)   : 0.7603

cat(sprintf("  No-skill baseline       : %.4f\n", prevalence))

##   No-skill baseline       : 0.2549

pr_l <- pr_curve(logit_proba, test_s$high_value); pr_l$model <- "Logistic"
pr_r <- pr_curve(rf_proba,    test$high_value);   pr_r$model <- "RF default"
pr_t <- pr_curve(tuned_proba, test$high_value);   pr_t$model <- "RF tuned"
pr_all <- rbind(pr_l, pr_r, pr_t)

ggplot(pr_all, aes(recall, precision, colour = model)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = prevalence, linetype = "dashed", colour = "grey50") +
  scale_colour_manual(values = c("Logistic" = "#2c7fb8",
                                 "RF default" = "#d95f0e",
                                 "RF tuned" = "#41b6c4")) +
  labs(title = "Precision-Recall Curve - High-Value class",
       x = "Recall", y = "Precision", colour = NULL) +
  theme_minimal() + theme(legend.position = c(0.8, 0.85))

Final Model Comparison

final <- data.frame(
  Model = c("Logistic Regression", "Random Forest (default)",
            "Random Forest (tuned)"),
  Accuracy  = round(c(m_logit$accuracy, m_rf$accuracy, m_tuned$accuracy), 3),
  Precision = round(c(m_logit$precision, m_rf$precision, m_tuned$precision), 3),
  `Recall (High)` = round(c(m_logit$recall, m_rf$recall, m_tuned$recall), 3),
  F1        = round(c(m_logit$f1, m_rf$f1, m_tuned$f1), 3),
  `ROC-AUC` = round(c(logit_auc, rf_auc, NA), 3),
  `PR-AUC`  = round(c(ap_logit, ap_rf, ap_tuned), 3),
  check.names = FALSE
)
knitr::kable(final, caption = "Final classification model comparison")

Final classification model comparison
Model	Accuracy	Precision	Recall (High)	F1	ROC-AUC	PR-AUC
Logistic Regression	0.818	0.608	0.800	0.691	0.892	0.738
Random Forest (default)	0.834	0.635	0.817	0.715	0.905	0.754
Random Forest (tuned)	0.834	0.636	0.816	0.715	NA	0.760

Conclusions for the report:

Results are checked under 5-fold cross-validation, so the reported F1 scores are stable rather than a lucky split.
Tuning the Random Forest gives the strongest overall trade-off, recovering high recall while keeping the forest’s better ranking.
On both ROC-AUC and PR-AUC the Random Forest beats logistic regression; PR-AUC is the fairer view here because the data is imbalanced.
The lazy baseline and the leakage proof earlier show why we report F1 / recall / AUC and why the leakage variables were excluded — the two core pitfalls of this task.

9. Discussion

The EDA shows that sale price is strongly right-skewed and that property size, borough, and building class are important factors. Gross square feet has the strongest numeric relationship with sale price, while property age contributes a weaker and non-linear pattern. These findings support the use of a log-transformed regression target and tree-based models.

For regression, Random Forest performs best, with RMSE 883,661.8, MAE 373,324.4, and R-squared 0.569. Decision Tree Regression performs better than Linear Regression, while the negative R-squared of the linear model indicates that a simple linear relationship is not sufficient for this dataset. Random Forest variable importance also identifies gross square feet, borough, and building class as major predictors.

For classification, both Random Forest models outperform Logistic Regression on the main comparison metrics. The default Random Forest achieves accuracy 0.834, F1 0.715, ROC-AUC 0.905, and PR-AUC 0.754. The tuned Random Forest keeps accuracy and F1 at 0.834 and 0.715 and produces the highest reported PR-AUC of 0.760. The cross-validation results also show a more stable F1 score for Random Forest than for Logistic Regression.

Location and property size are important because property values differ substantially across boroughs and because larger buildings generally sell for higher prices. Tree-based models perform better because they can represent non-linear relationships and interactions, such as different size effects across locations and building classes.

The project has several limitations. The data covers only 2016-2017 and does not include all factors that may influence property value, such as market conditions, renovation status, property quality, exact street-level location, and neighborhood amenities. The High Value classification target is also based on a dataset-specific price threshold, so its interpretation may change in another period or market.

10. Conclusion

This project followed a complete data science workflow, including data cleaning, feature engineering, exploratory analysis, regression modelling, and classification modelling.

The first modelling question can be answered positively: NYC property sale prices can be predicted to some extent using property characteristics. Random Forest Regression is the best regression model in the existing results because it has the lowest RMSE and MAE and the highest R-squared among the three models.

The second modelling question can also be answered positively: high-value property transactions can be classified using property and transaction characteristics. Random Forest provides the strongest classification performance, with the tuned model giving the best reported overall precision-recall trade-off.

Future work could include more detailed location information, market indicators, renovation and property-quality variables, and data from additional years. These additions may improve prediction accuracy and make the models more useful across different market conditions.

11. References

Kaggle. NYC Property Sales Dataset. https://www.kaggle.com/datasets/new-york-city/nyc-property-sales
New York City Department of Finance. Rolling Sales Data.
R Core Team. R: A Language and Environment for Statistical Computing.
R packages used in the project include tidyverse, readr, ggplot2, scales, corrplot, patchwork, rpart, randomForest, knitr, and rmarkdown.

Data-Driven Insights into NYC Property Market: Price Prediction and Transaction Classification

WQD7004 Group Project

June 2026