<jupyter_start><jupyter_text>**Linear Regression to Predict Housing Sale Price** **Deep Dive Data Science Project 2**

---
<NAME> 3/3/2021  Problem Definition We desire to predict home sale prices based on various features of the property.

Output to be produced: Predicted sale price, in US dollars.

Additional goal: Attempt to determine a minimal set of features.  Supporting Functions<jupyter_code>def root_mean_squared_percentage_error(y_true, y_pred):
    """Signature(y_true, y_pred). Returns float representing RMSE in percent."""
    rmspe = np.sqrt(np.mean(np.square((y_true-y_pred)/y_true)))*100
    return rmspe<jupyter_output><empty_output><jupyter_text>Data Collection<jupyter_code># Import required packages

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import sklearn.model_selection as model_selection
from sklearn import metrics

from IPython.display import Image

from sklearn.linear_model import Ridge
# The data were sourced from the instructor.  It contains information related to home sale transactions.
# File Path Settings

file_name = '/content/drive/MyDrive/Data Science Student Material/Project 2 - Cohort 4/Data/Housing Data.csv'

train = pd.read_csv(file_name)

train.info()
print(train.head())
print(train.shape)<jupyter_output>PID  MS SubClass MS Zoning  ...  Sale Type  Sale Condition SalePrice
0  905101070           20        RL  ...        WD           Normal    115400
1  905101330           90        RL  ...        WD           Normal     90000
2  903454090           50        RM  ...        WD           Normal    141000
3  533244030           60        FV  ...        WD           Normal    254000
4  909252020           70        RL  ...        WD           Normal    155000

[5 rows x 81 columns]
(2637, 81)<jupyter_text>Data Cleaning<jupyter_code># Drop spurious columns

train_clean = train.copy()

train_clean.drop('PID', axis=1, inplace=True)

train_clean.head()
# MS_SubClass is coded as an integer, but really represents different categories,
# so change to String type

train_clean['MS SubClass'] = train_clean['MS SubClass'].apply(str)

#train_clean.info()
# Replace 'Lot Frontage' NAs with zero (0.0)

train_clean['Lot Frontage'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Alley' NAs with 'No_Alley'

train_clean['Alley'].fillna(value='No_Alley', inplace=True)

#train_clean.info()
# Replace 'Mas Vnr Type' NAs with 'No_Vnr'

train_clean['Mas Vnr Type'].fillna(value='No_Vnr', inplace=True)

#train_clean.info()
# Replace 'Mas Vnr Area' NAs with 0.0

train_clean['Mas Vnr Area'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Bsmt Qual' NAs with 'No_Bsmt'

train_clean['Bsmt Qual'].fillna(value='No_Bsmt', inplace=True)

#train_clean.info()
# Replace 'Bsmt Cond' NAs with 'No_Bsmt'

train_clean['Bsmt Cond'].fillna(value='No_Bsmt', inplace=True)

#train_clean.info()
# Replace 'Bsmt Exposure' NAs with 'No_Bsmt'

train_clean['Bsmt Exposure'].fillna(value='No_Bsmt', inplace=True)

#train_clean.info()
# Replace 'BsmtFin Type 1' NAs with 'No_Bsmt'

train_clean['BsmtFin Type 1'].fillna(value='No_Bsmt', inplace=True)

#train_clean.info()
# Replace 'BsmtFin SF 1' NAs with 0.0

train_clean['BsmtFin SF 1'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'BsmtFin Type 2' NAs with 'No_Bsmt'

train_clean['BsmtFin Type 2'].fillna(value='No_Bsmt', inplace=True)

#train_clean.info()
# Replace 'BsmtFin SF 2' NAs with 0.0

train_clean['BsmtFin SF 2'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Bsmt Unf SF' NAs with 0.0

train_clean['Bsmt Unf SF'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Total Bsmt SF' NAs with 0.0

train_clean['Total Bsmt SF'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Bsmt Full Bath' NAs with 0.0

train_clean['Bsmt Full Bath'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Bsmt Half Bath' NAs with 0.0

train_clean['Bsmt Half Bath'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Fireplace Qu' NAs with 'No_FirePl'

train_clean['Fireplace Qu'].fillna(value='No_FirePl', inplace=True)

#train_clean.info()
# Replace 'Garage Type' NAs with 'No_Garage'

train_clean['Garage Type'].fillna(value='No_Garage', inplace=True)

#train_clean.info()
# Decision made to drop 'Garage Yr Blt' column
# due to low correlation and no intuitive NA replacement.

train_clean.drop('Garage Yr Blt', axis=1, inplace=True)

train_clean.head()
# Replace 'Garage Finish' NAs with 'No_Garage'

train_clean['Garage Finish'].fillna(value='No_Garage', inplace=True)

#train_clean.info()
# Replace 'Garage Cars' NAs with 0.0

train_clean['Garage Cars'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Garage Area' NAs with 0.0

train_clean['Garage Area'].fillna(value=0.0, inplace=True)

#train_clean.info()
# Replace 'Garage Qual' NAs with 'No_Garage'

train_clean['Garage Qual'].fillna(value='No_Garage', inplace=True)

#train_clean.info()
# Replace 'Garage Cond' NAs with 'No_Garage'

train_clean['Garage Cond'].fillna(value='No_Garage', inplace=True)

#train_clean.info()
# Replace 'Pool QC' NAs with 'No_Pool'

train_clean['Pool QC'].fillna(value='No_Pool', inplace=True)

#train_clean.info()
# Replace 'Fence' NAs with 'No_Fence'

train_clean['Fence'].fillna(value='No_Fence', inplace=True)

#train_clean.info()
# Replace 'Misc Feature' NAs with 'No_Misc'

train_clean['Misc Feature'].fillna(value='No_Misc', inplace=True)

#train_clean.info()
# Removal of non-'Normal' 'Sale Condition' observations after EDA

train_clean = train_clean[train_clean['Sale Condition'] == 'Normal']

train_clean.shape  # Correct number of rows; 2637 - 471
# Drop misc observations: 2001 - Basement incomplete; 162 - Garage Yr Blt inconsistent; 1555 - only Tennis Ct

train_clean = train_clean.drop([2001, 162, 1555], axis=0)

train_clean.shape
train_clean = train_clean.reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis<jupyter_code># Examine SalePrice histogram

plt.hist(train_clean['SalePrice'], bins=40);
plt.title("Sale Price Histogram", size=20);
plt.xlabel("Actual Sale Price, US$", size=15)
plt.ylabel("Number of Sales", size=15);

# Separate out numeric, non-numeric/object features

train_numeric = train_clean.select_dtypes(include=['int64', 'float64'])

train_object = train_clean.select_dtypes(include='object')
numeric_feature_corr = train_numeric.corr()['SalePrice'].sort_values(ascending=False)[1:]

print(numeric_feature_corr.head())
print(numeric_feature_corr.tail(), '\n')

golden_features = numeric_feature_corr.abs().sort_values(ascending=False)
golden_features
import statsmodels.api as sm

actual_outcomes = train_numeric['SalePrice'].copy()

train_numeric_one_feature = pd.DataFrame(train_numeric[golden_features.index[0]].copy())

train_numeric_one_feature['Intercept'] = 1

train_numeric_one_feature
# Instantiate, fit, examine model results

model_one_feature = sm.OLS(actual_outcomes, train_numeric_one_feature).fit()

prediction_one_feature = model_one_feature.predict(train_numeric_one_feature)

print(model_one_feature.summary())

plt.figure(figsize=(12,10))
plt.scatter(actual_outcomes, prediction_one_feature);
plt.plot([0, max(prediction_one_feature)], [0, max(prediction_one_feature)], c='red');
# Try two features

train_numeric_two_features = pd.DataFrame(train_numeric[golden_features.index[:2]].copy())
train_numeric_two_features['Intercept'] = 1
train_numeric_two_features
model_two_features = sm.OLS(actual_outcomes, train_numeric_two_features).fit()
prediction_two_features = model_two_features.predict(train_numeric_two_features)
print(model_two_features.summary())
plt.figure(figsize=(12,10))
plt.scatter(actual_outcomes, prediction_two_features);
plt.plot([0, max(prediction_two_features)], [0, max(prediction_two_features)], c='red');
# We see some outliers, so we investigate

prediction_two_features[prediction_two_features>450000]
#train_clean.iloc[1791,:]  # Sale Condition == Partial
#train_clean.iloc[1939,:]  # Sale Condition == Partial
#train_clean.iloc[2494,:]  # Sale Condition == Abnormal
train_clean[train_clean['Sale Condition'] != 'Normal']  # 471 observations; 17.9% of data

# Removal of non-Normal observations performed above

model_all_features = sm.OLS(actual_outcomes, train_numeric[golden_features.index[:33]]).fit()
prediction_all_features = model_all_features.predict(train_numeric[golden_features.index[:33]])
print(model_all_features.summary())
plt.figure(figsize=(10, 10))
plt.scatter(actual_outcomes, prediction_all_features);
plt.plot([0, max(prediction_all_features)], [0, max(prediction_all_features)], c='red');
root_mean_squared_percentage_error(actual_outcomes, prediction_all_features)<jupyter_output><empty_output><jupyter_text>Data Processing  Cross Validation<jupyter_code># Run a split/fit to assess initial accuracy

X_train, X_test, y_train, y_test = model_selection.train_test_split(train_numeric[golden_features.index[:33]], actual_outcomes, test_size=0.20)

model_all_features_split = sm.OLS(y_train, X_train).fit()

prediction_all_features_split = model_all_features_split.predict(X_test)

root_mean_squared_percentage_error(y_test, prediction_all_features_split)
# Run cross-validation over many iterations to assess average RMS%, stddev

numLoops = 100
results  = np.zeros(numLoops)
resArray = np.zeros((2,2))

for idx in range(numLoops):
  X_train, X_test, y_train, y_test = model_selection.train_test_split(train_numeric[golden_features.index[:33]], actual_outcomes, test_size=0.20)
  model_all_features_split = sm.OLS(y_train, X_train).fit()
  prediction_all_features_split = model_all_features_split.predict(X_test)

  results[idx] =root_mean_squared_percentage_error(y_test, prediction_all_features_split)

print(results.mean())
print(results.std())<jupyter_output>15.293904057168426
1.303298814440501<jupyter_text>Scale Numerical Features<jupyter_code>from sklearn import preprocessing

scaler = preprocessing.MinMaxScaler()  # RobustScaler underperformed

train_numeric_scaled = scaler.fit_transform(train_numeric)

train_numeric_scaled = pd.DataFrame(train_numeric_scaled, columns=train_numeric.columns)
train_numeric_scaled.describe()
train_numeric_scaled.shape
train_numeric_scaled<jupyter_output><empty_output><jupyter_text>Categorical Features<jupyter_code>train_object['Neighborhood']
# One-hot encode categorical features

train_one_hot = pd.get_dummies(train_object['Neighborhood'])

train_object['Neighborhood'].value_counts()
print(train_one_hot.shape)
train_one_hot
train_numeric_one_hot = pd.concat([train_numeric_scaled, train_one_hot], axis=1)
train_numeric_one_hot.drop(['SalePrice'], axis=1, inplace=True)
train_numeric_one_hot['Intercept'] = 1

train_numeric_one_hot
model_one_hot = sm.OLS(actual_outcomes, train_numeric_one_hot).fit()

predict_one_hot = model_one_hot.predict(train_numeric_one_hot)

root_mean_squared_percentage_error(actual_outcomes, predict_one_hot)
# Loop over all Categorical features to see most impactful

results = pd.DataFrame(index=['RMS%'])
train_numeric_scaled.drop(['SalePrice'], axis=1, inplace=True)
train_numeric_scaled['Intercept'] = 1

for feature in train_object.columns:
  train_one_hot = pd.get_dummies(train_object[feature])
  train_numeric_one_hot = pd.concat([train_numeric_scaled, train_one_hot], axis=1)
  model_one_hot = sm.OLS(actual_outcomes, train_numeric_one_hot).fit()
  predict_one_hot = model_one_hot.predict(train_numeric_one_hot)
  results[feature] = root_mean_squared_percentage_error(actual_outcomes, predict_one_hot)

results = results.sort_values(by='RMS%', axis=1)
results<jupyter_output><empty_output><jupyter_text>Add top categorical features<jupyter_code>train_one_hot = pd.get_dummies(train_object[results.columns[:7]])
train_numeric_one_hot = pd.concat([train_numeric_scaled, train_one_hot], axis=1)
model_one_hot = sm.OLS(actual_outcomes, train_numeric_one_hot).fit()
predict_one_hot = model_one_hot.predict(train_numeric_one_hot)
root_mean_squared_percentage_error(actual_outcomes, predict_one_hot)
# Run cross-validation over many iterations to assess average RMS%, stddev

numLoops = 100
results  = np.zeros(numLoops)
resArray = np.zeros((2,2))

for idx in range(numLoops):
  X_train, X_test, y_train, y_test = model_selection.train_test_split(train_numeric_one_hot, actual_outcomes, test_size=0.20)
  model_all_features_split = sm.OLS(y_train, X_train).fit()
  prediction_all_features_split = model_all_features_split.predict(X_test)

  results[idx] =root_mean_squared_percentage_error(y_test, prediction_all_features_split)

print(results.mean())
print(results.std())<jupyter_output>11.165483043264178
0.716631590709585<jupyter_text>Non-linear Transform Attempt to make predictions more accurate by transforming outcomes via nonlinear (quadratic) function.<jupyter_code>degree = 0.3
nonlinear_outcomes = actual_outcomes.pow(degree)

model_one_hot = sm.OLS(nonlinear_outcomes, train_numeric_one_hot).fit()
predict_one_hot = model_one_hot.predict(train_numeric_one_hot)

real_predicted_outcomes = predict_one_hot.pow(1/degree)

root_mean_squared_percentage_error(actual_outcomes, real_predicted_outcomes)

plt.figure(figsize=(12,10));
plt.scatter(actual_outcomes, real_predicted_outcomes);
plt.plot([0, max(real_predicted_outcomes)],[0, max(real_predicted_outcomes)], c='red');<jupyter_output><empty_output><jupyter_text>Determine optimal degree<jupyter_code>degrees = np.arange(.1, 1, .01)

percentage_error_loop = np.zeros(len(degrees))

for idx, degree in enumerate(degrees):
  
  nonlinear_outcomes = actual_outcomes.pow(degree)

  model_loop = sm.OLS(nonlinear_outcomes, train_numeric_one_hot).fit()
  predict_loop = model_loop.predict(train_numeric_one_hot)

  real_predicted_outcomes_loop = predict_loop.pow(1/degree)

  percentage_error_loop[idx] = root_mean_squared_percentage_error(actual_outcomes, real_predicted_outcomes_loop)
plt.plot(percentage_error_loop)
percentage_error_loop.argmin()/100  # Optimal degree<jupyter_output><empty_output><jupyter_text>Incorporate optimal degree, showing ~ 1.1% improvement<jupyter_code>degree = 0.22
nonlinear_outcomes = actual_outcomes.pow(degree)

model_one_hot = sm.OLS(nonlinear_outcomes, train_numeric_one_hot).fit()
predict_one_hot = model_one_hot.predict(train_numeric_one_hot)

real_predicted_outcomes = predict_one_hot.pow(1/degree)

print(root_mean_squared_percentage_error(actual_outcomes, real_predicted_outcomes))

plt.figure(figsize=(12,10));
plt.scatter(actual_outcomes, real_predicted_outcomes);
plt.plot([0, max(real_predicted_outcomes)],[0, max(real_predicted_outcomes)], c='red');<jupyter_output>9.280380378615785<jupyter_text>Ridge Regression<jupyter_code>model_ridge = Ridge(alpha=0.001)
model_ridge.fit(train_numeric_one_hot, nonlinear_outcomes)
predict_ridge = model_ridge.predict(train_numeric_one_hot)
predict_ridge_df = pd.Series(predict_ridge)
real_predicted_ridge = predict_ridge_df.pow(1/degree)

root_mean_squared_percentage_error(actual_outcomes, real_predicted_ridge)<jupyter_output><empty_output><jupyter_text>Final feature dataframe<jupyter_code>train_all_one_hot = pd.get_dummies(train_object)
train_ridge = pd.concat([train_numeric, train_all_one_hot], axis=1)
train_ridge.drop(['SalePrice'], axis=1, inplace=True)
train_ridge.head()
# Find optimal alpha

alpha_value = [1e-3, 1e-2, 3e-2, 6e-2, 1e-1, 3e-1, 6e-1, 1, 10]

percentage_error_ridge = np.zeros(len(alpha_value))

for idx, alpha_val in enumerate(alpha_value):
  ridgereg = Ridge(alpha=alpha_val)
  ridgereg.fit(train_ridge, nonlinear_outcomes)
  pred_ridge = pd.Series(ridgereg.predict(train_ridge))
  percentage_error_ridge[idx] = root_mean_squared_percentage_error(actual_outcomes, pred_ridge.pow(1/degree))

percentage_error_ridge
# Select best alpha

best_alpha = 0.1

plt.plot(percentage_error_ridge)
# Cross-validation

numLoops = 100
results  = np.zeros(numLoops)

for idx in range(numLoops):
  X_train, X_test, y_train, y_test = model_selection.train_test_split(train_ridge, nonlinear_outcomes, test_size=0.20)
  model_ridge = Ridge(alpha = best_alpha)
  model_ridge.fit(X_train, y_train)
  pred_ridge = pd.Series(model_ridge.predict(X_test))
  results[idx] = root_mean_squared_percentage_error(y_test, pred_ridge)

print(results.mean())
print(results.std())
model_ridge = Ridge(alpha=0.001)
model_ridge.fit(train_numeric_one_hot, nonlinear_outcomes)
predict_ridge = model_ridge.predict(train_numeric_one_hot)
predict_ridge_df = pd.Series(predict_ridge)
real_predicted_ridge = predict_ridge_df.pow(1/degree)
plt.figure(figsize=(8, 6))
plt.scatter(actual_outcomes, real_predicted_ridge);
plt.plot([0, max(actual_outcomes)], [0, max(actual_outcomes)], c='red');<jupyter_output><empty_output><jupyter_text>Data Visualizations  Predicting Blind Test Data<jupyter_code>best_degree = 0.22
best_alpha = 0.1

final_actual_outcomes = train_clean['SalePrice'].copy()
final_nonlinear_outcomes = final_actual_outcomes.pow(degree)

final_numeric = train_numeric.drop(['SalePrice'], axis=1).copy()

final_category = train_object.copy()

final_all_features = pd.concat([final_numeric, final_category])
final_one_hot = pd.get_dummies(final_category)
final_all_one_hot = pd.concat([final_numeric, final_one_hot], axis=1)
final_all_one_hot.shape<jupyter_output><empty_output><jupyter_text>Expected Prediction Accuracy<jupyter_code># Cross-validation

numLoops = 100
results  = np.zeros(numLoops)

for idx in range(numLoops):
  X_train, X_test, y_train, y_test = model_selection.train_test_split(final_all_one_hot, final_nonlinear_outcomes, test_size=0.20)
  model_ridge = Ridge(alpha = best_alpha)
  model_ridge.fit(X_train, y_train)
  pred_ridge = pd.Series(model_ridge.predict(X_test))
  results[idx] = root_mean_squared_percentage_error(y_test, pred_ridge)

print(results.mean())
print(results.std())
plt.hist(results, bins=10)<jupyter_output><empty_output><jupyter_text>Train final model with all data<jupyter_code>model_final = Ridge(alpha=best_alpha)
model_final.fit(final_all_one_hot, final_nonlinear_outcomes)<jupyter_output><empty_output><jupyter_text>Load test data<jupyter_code>test_data = pd.read_csv('/content/drive/MyDrive/Data Science Student Material/Project 2 - Cohort 4/Data/Housing Data Blind Test.csv')
test_data.shape<jupyter_output><empty_output><jupyter_text>Adjust test data to match training set<jupyter_code># Cast MS SubClass as String
test_data['MS SubClass'] = test_data['MS SubClass'].apply(str)

# Check Sale Condition for Normal label
print(test_data['Sale Condition'].value_counts())

# Check one_hot_encoding labels
test_data_numeric = test_data[final_numeric.columns]
test_data_category = test_data[final_category.columns]
test_category_one_hot = pd.get_dummies(test_data_category)
print('Examine test_category_one_hot.shape vs. final_one_hot.shape:', test_category_one_hot.shape, final_one_hot.shape)<jupyter_output>Normal    246
Name: Sale Condition, dtype: int64
Examine test_category_one_hot.shape vs. final_one_hot.shape: (246, 228) (2163, 282)<jupyter_text>We see there are missing labels (one_hot categories) in the test data, so adjust for the discrepancy.<jupyter_code>test_data_one_hot = pd.concat([test_data_numeric, test_category_one_hot], axis=1)

# Create new df with correct size/indexing
test_data_full = pd.DataFrame(0, index=np.arange(test_data.shape[0]), columns=final_all_one_hot.columns)
test_data_full.head()
# Form intersection of available test-one-hot columns with full data one-hot columns
column_intersect = final_all_one_hot.columns.intersection(test_data_one_hot.columns)
column_intersect.size
# Fill new df with available test data/one-hot labels
test_data_full[column_intersect] = test_data_one_hot[column_intersect]<jupyter_output><empty_output><jupyter_text>Predict test dataset<jupyter_code>predict_blind_test = pd.Series(model_final.predict(test_data_full))
predict_final = predict_blind_test.pow(1/best_degree)
predict_final
# Write final predicted data to final_submission.csv

#predict_final.to_csv('/content/drive/MyDrive/Exercises - Ken/Project 2/final_submission.csv', index=False)<jupyter_output><empty_output><jupyter_text>Conclusions - Data supports Regression approach
- Accuracy improvements with additional techniques
- Cross-validated accuracy 88.5%
- Model stable<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10908212560386474, 'fraction_numerical': 0.02323671497584541, 'mean_word_length': 6.489507959479016, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '558213', 'license': 'MIT', 'n_tokens_mistral': 7848, 'n_tokens_neox': 7236, 'n_words': 1294, 'path': 'Project_2.ipynb', 'repo_name': 'kengustafson/DeepDive-DataScience', 'size': '20707'}