<jupyter_start><jupyter_text>Part B: Zillow House Price Regression  Variable Description:- Price - price of property (dependent variable)- Area - size of the property (square feet)- Bedrooms - number of bedrooms in the property- Bathrooms - number of bathrooms in the property- Stories - number of floors in the property- Mainroad - indicates if the property is located close to a main road- Guestroom - indicates if the property has a guest room- Basement - indicates if the property has a basement- Hot Water Heating - indicates if the property uses gas for hot water heater- Air Conditioning - indicates if the property has an air conditioning unit- Parking - number of parking spaces- Pref Area - indicates if the property is located in a preferred neighbourhood of the city- Furninghing - indicates if the house is furnished or not   Load and Explore Data Structure<jupyter_code># Import necessary packages
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# Load the dataset
raw_zillow = pd.read_csv("House_Price.csv")
raw_zillow.head()
# Rename columns for easier coding
raw_zillow = raw_zillow.rename(columns = {"bathrooms":"bath", "hotwaterheating":"hotwater", "airconditioning": "aircon",
                      "furnishingstatus":"furnishing"})
raw_zillow.columns
# View structure of the dataset
raw_zillow.info()
# View summary statistics
raw_zillow.describe()
# Check for any null values
raw_zillow.isnull().sum()<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis<jupyter_code># Distribution of Price (dependent variable)
sns.boxplot(x=raw_zillow.price, data=raw_zillow);
plt.title(" Boxplot of Price - Slightly right skewed", fontsize=13)
plt.ticklabel_format(style='plain', axis='x')
sns.set_style("darkgrid")
plt.figure(figsize=(18,5))
plt.subplots_adjust(wspace = 0.5)


plt.subplot(131)
sns.countplot(x="bedrooms", data=raw_zillow)
plt.ticklabel_format(style='plain', axis='y')
plt.title("Count of Bedrooms", fontsize=13)


plt.subplot(132)
sns.barplot(x="bedrooms", y="price", data=raw_zillow, errwidth=False)
plt.ticklabel_format(style='plain', axis='y')
plt.title("Average Price by Bedrooms", fontsize=13);

plt.subplot(133)
sns.boxplot(x="bedrooms", y="price", data=raw_zillow)
plt.ticklabel_format(style='plain', axis='y')
plt.title("Boxplot of Price by Bedrooms", fontsize=13);<jupyter_output><empty_output><jupyter_text>* Observation - there are alot more 3 bedroom properties. Average price of 3, 4, and 5 bedrooms(the other variables are also involved in the price) are much higher than others. Additionally 3, 4, and 5 bedroom properties have a wider range of prices. * The charts below examine the other categorical variables by show the frequency along with the average price.<jupyter_code># Function to compare frequency and average price of categorical variables
def my_function(df, var1, var2):
    
    sns.set_style("darkgrid")
    plt.figure(figsize=(18,5))
    plt.subplots_adjust(wspace = 0.5)
   
    title = str(var1).title()
    
    plt.subplot(131)
    sns.countplot(x=str(var1), data=df)
    plt.ticklabel_format(style='plain', axis='y')
    plt.title("Count ({})".format(title), fontsize=13)
    
    plt.subplot(132)
    sns.barplot(x=str(var1), y=str(var2), data=df, errwidth=False)
    plt.ticklabel_format(style='plain', axis='y')
    plt.title("Average Price ({})".format(title), fontsize=13)
    
    plt.subplot(133)
    sns.boxplot(x=str(var1), y=str(var2), data=df)
    plt.ticklabel_format(style='plain', axis='y')
    plt.title("Price Distribution ({})".format(title), fontsize=13);
my_function(raw_zillow, "mainroad", "price");
my_function(raw_zillow, "basement", "price");
my_function(raw_zillow, "hotwater", "price");
my_function(raw_zillow, "guestroom", "price");
my_function(raw_zillow, "hotwater", "price");
my_function(raw_zillow, "guestroom", "price");
my_function(raw_zillow, "aircon", "price");
my_function(raw_zillow, "parking", "price");
my_function(raw_zillow, "prefarea", "price");
my_function(raw_zillow, "furnishing", "price");<jupyter_output><empty_output><jupyter_text>transform "yes" and "no" values to 1 and 0<jupyter_code># iterate over the columns
for col in raw_zillow.columns:
    if raw_zillow.loc[0, col] in ["yes", "no"]: # if the first row of each column corresponds to yes or no
        raw_zillow[col] = raw_zillow[col].map({"yes":1, "no":0}) # convert them to 1 or 0 

raw_zillow.head(10)<jupyter_output><empty_output><jupyter_text>convert the categorical variable "furnishingstatus" to a dummy variable<jupyter_code># convert it to a dummy variable and drop the first dummy so  [semi and unfurnished are 0] = "furnished" 
raw_zillow = pd.get_dummies(raw_zillow, columns = ["furnishing"], drop_first = True)
raw_zillow.head(10)<jupyter_output><empty_output><jupyter_text>Split the dataset into train and test sets<jupyter_code>from sklearn.model_selection import train_test_split

y = raw_zillow["price"] # dependent variable = "price"
X = raw_zillow.iloc[:, 1:] # independent variables = the rest of the variables

x_train, x_test, y_train, y_test = train_test_split(X, y, train_size = 0.8, random_state = 123) # train:test = 8:2
import statsmodels.api as sm

x_train_const = sm.add_constant(x_train) # add a constant column

train_ml_result = sm.OLS(y_train, x_train_const).fit()
train_ml_result.summary()<jupyter_output><empty_output><jupyter_text>* F-test p-value of the multivariate regression model presents the statistical significance of the overall model and an approximate 68% of the variation in "Yearly-Amount Spent" in population is explained by the predictors based on the adjusted R-squred.  * Each of the independent variables excepting "bedrooms" & "semi-furnished" seems to have their unique ability to predict the "housing price" based on their t-test p-value smaller than 0.05.  * Create a annotated heatmap to check correlations among variables<jupyter_code># concatenate x_train and y_train to make a single dataframe

train_df = pd.concat([y_train, x_train], axis = 1)
train_df.head()
import matplotlib.pyplot as plt
import seaborn as sns

fig = plt.figure(figsize = (14, 10))
sns.heatmap(train_df.corr(), annot = True)<jupyter_output><empty_output><jupyter_text>* no multicollinearity issue and no variables significantly correlated with "bedrooms" or "semi-furnished" are observed. That is, droppingn both or either the variables will unlikely affect other variables * Determine whether to drop the "bedrooms" & "semi-furnished" variables<jupyter_code># First, explore the "bedrooms" variable

x_train["bedrooms"].describe()
# histogram for "bedrooms"
plt.hist(x_train["bedrooms"])<jupyter_output><empty_output><jupyter_text>* mean and median is highly close and the histogram also shows the "bedrooms" is not skewed * compare the adj r-squared before and after the variable is dropped from the dataset<jupyter_code># remove the "bedrooms" variable from x_train 
x_train2 = x_train.drop(["bedrooms"], axis = 1)
x_train2.head()
# fit the mlr model again with the free-bedrooms train dataset

train_ml_result2 = sm.OLS(y_train, sm.add_constant(x_train2)).fit()
train_ml_result2.summary()<jupyter_output><empty_output><jupyter_text>* As expected, there's no change in adjusted r-squared * Drop the semi-furnished<jupyter_code>x_train2 = x_train2.drop(["furnishing_semi-furnished"], axis = 1)
x_train2.head()
# fit the mlr model again with the new train dataset

train_ml_result3 = sm.OLS(y_train, sm.add_constant(x_train2)).fit()
train_ml_result3.summary()<jupyter_output><empty_output><jupyter_text>* after dropping the "semi-furnished" variable, the adj r-squared has increased by 0.001 and all t-test p-values and F-test p-value are smaller than 0.05, which indicate that we can reasonbly be sure that the slope of each predictor is not 0 in population holding the other predictors constant and the overall model is also statistically significant  Diagnose the model<jupyter_code># Residual histogram
train_resid = train_ml_result3.resid # get the residuals of the MLR model
plt.hist(train_resid, bins=20, rwidth=0.9)
plt.title("Residual Histogram")
plt.xlabel("Residual")<jupyter_output><empty_output><jupyter_text>* the residual histogram seems to have ourliers  * detect and remove outliers  1) conduct a normalized residual test to detect outliers<jupyter_code># whether the normalized residuals are lying within a rule-of-thumb threshold 
norm_resid = train_ml_result3.get_influence().resid_studentized_internal

plt.figure(figsize = (14, 4))
plt.stem(norm_resid)
plt.axhline(2, c="g", ls="--") # outliers if the normalized residuals are bigger than |2| ~ |4|
plt.axhline(-2, c="g", ls="--") # set the threshold as |2|
plt.xlabel("Normalized Residuals")
# Detect the indices of norm_resid values outbound the range (|3|)
import numpy as np

norm_resid_df = pd.DataFrame(norm_resid, columns = ["norm_resid"])

resid_outlier = norm_resid_df[np.abs(norm_resid_df["norm_resid"]) > 2].index # get the indices of the observations lying outside the range
# create a data frame with x_train & y_train and reset the index to match the detected outlier indices 

train_df2 = pd.concat([y_train, x_train2], axis = 1)
train_df2 = train_df2.reset_index(drop=True)
train_df2.head()
# drop the detected indexes
train_df3 = train_df2.drop(resid_outlier)
train_df3.head()<jupyter_output><empty_output><jupyter_text>* Conduct a z-score test too with the new data frame created right above<jupyter_code>from scipy import stats

# remove values whose z-score > |3| across all the columns 
train_df4 = train_df3[(np.abs(stats.zscore(train_df3)) < 3).all(axis=1)] # extract the only rows whose absolute z-scores are smaller than 3
train_df4 = train_df4.reset_index(drop = True)
train_df4.head()
# extract the y and x variables from the train_df4 
# and then fit a mlr model again with the new datasets
y_train = train_df4["price"]
x_train = train_df4.iloc[:, 1:]
train_ml_result4 = sm.OLS(y_train, sm.add_constant(x_train)).fit()<jupyter_output><empty_output><jupyter_text>* Check the Cook's distance<jupyter_code># create a graph for cook's distance
train_influence = train_ml_result4.get_influence()
(c, p) = train_influence.cooks_distance
plt.stem(np.arange(len(x_train)),
np.round(c, 3),
markerfmt=",")
plt.title("Cooks distance")
plt.xlabel("Row index")
plt.ylabel("Cooks Distance")
train_cooks = train_influence.cooks_distance[0] # get cook's distance
outlier_list = [] 
for i in range(len(train_cooks)):
    if train_cooks[i] > 4/(len(train_cooks)-len(train_ml_result4.params)-1): # use the rule of thumb suggested by Fox to find large D
        outlier_list.append(i)  # put the filtered indices into the empty list
print(outlier_list) 
train_df5 = train_df4.drop(outlier_list) # remove the rows matched to the indices in the list
train_df5 = train_df5.reset_index(drop=True) # reset the index
y_train = train_df5["price"] # extract the dependent variable from the data frame
x_train = train_df5.iloc[:, 1:] # independent variables 
train_ml_result4 = sm.OLS(y_train, sm.add_constant(x_train)).fit() # re-fit<jupyter_output><empty_output><jupyter_text>* Create a scatter plot with the predicted y values and studentized residuals<jupyter_code># Scatter plot with studentized residuals and predicted y values
norm_resid2 = train_ml_result4.get_influence().resid_studentized_internal # normalized the residuals
y_predict = train_ml_result4.fittedvalues # get predicted y values

fig = plt.figure(figsize=(14, 9))
ax = fig.add_subplot()
ax.scatter(y_predict, norm_resid2, alpha = 0.6)
sns.regplot(y_predict, norm_resid2, scatter = False, ci = False, lowess = True,  # add a line showing the shape of the scatter plot
          line_kws={"color":"red", "lw":1.2})
ax.axhline(y=0, c="black", lw=1.3, ls = "--")
ax.set_title("Studentized Residual Vs Predicted y")
ax.set_xlabel("Predicted y")
ax.set_ylabel("Studentized Residual")<jupyter_output><empty_output><jupyter_text>* neither serious heteroscedasticity nor non-normality features are observed.  * create a Q-Q plot to make sure the model's normality<jupyter_code>sm.qqplot(norm_resid2, line = "45") 
plt.title("Q-Q plot")
plt.show()<jupyter_output><empty_output><jupyter_text>* the data is well matched to the normal distribution values of the corresponding quantiles  * then, standardize the non-categorical independent variables to compare their impacts on the Housing Price<jupyter_code>from sklearn.preprocessing import StandardScaler  # load it to standardize independent variables 
scaler = StandardScaler()

part1 = x_train.iloc[:, [0, 1, 2, 8]] # extract only non-categorical columns
part2 = x_train.iloc[:, [3, 4, 5, 6, 7, 9, 10]] # part2 will be concatenated with part1 after scaling

part1_df = pd.DataFrame(scaler.fit_transform(part1), columns = part1.columns) # scaling part1

scaled_x_train = pd.concat([part1_df, part2], axis = 1) # combine the standardaized columns and part2 into a single data frame
scaled_x_train.head()
scaled_result = sm.OLS(y_train, sm.add_constant(scaled_x_train)).fit()
scaled_result.summary()
# resulting coefficients of unscaled predictors
from sklearn.linear_model import LinearRegression 
train_result_sk = LinearRegression().fit(x_train, y_train)
train_result_sk.coef_<jupyter_output><empty_output><jupyter_text>* after removing outliers, the model has been improved : 1) Based on the adjusted r-squared, the strength of the overall relationship with dependent variable has increased from 0.681 to      0.75. Now, approximately 79.3% of the variations in "Price" is predictable by the collective independent variables in population.  2)  the p-value of the intercept has decreased below 0.05 3) F-test P-value has also decreased a huge amount -> on the basis of the coefficients of each predictor, "prefarea" shows the stongest positive impacts on the housing price with the interpretation as the house at a prefering area increasing the price by about 689,135 dollars -> the top 5 indipendent variables can be listed up by their affects on the dependent variable (price): (1) prefarea (2) aircon (3) area (4) stories (5) furnishing status - furnished * in the outlier-free model, "hotwaterheating" has become statistically insignificant individually.  * drop the variable from the model and yield the finall results<jupyter_code>x_train = x_train.drop(["hotwater"], axis = 1)
final_result_sk = LinearRegression().fit(x_train, y_train)
final_result = sm.OLS(y_train, sm.add_constant(x_train)).fit()
final_result.summary()<jupyter_output><empty_output><jupyter_text>Test the trained model with the test data<jupyter_code>x_test = x_test.drop(["bedrooms", "hotwater", "furnishing_semi-furnished"], axis = 1) # drop the three columns from x_test
from sklearn.metrics import r2_score, mean_squared_error
pred_y_test = final_result_sk.predict(x_test) # get the predicted y values by plugging x_test values 
np.abs(r2_score(y_test, pred_y_test)) # get the r_squred score
# scatter plot with the predicted y and the test y 
plt.scatter(y_test, pred_y_test)
plt.title("Predicted y Vs Test y")
plt.xlabel("Test y")
plt.ylabel("Predicted y")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10386616610683001, 'fraction_numerical': 0.015609563327405651, 'mean_word_length': 5.971533516988062, 'pattern_counts': {'":': 8, '<': 78, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '150393', 'license': 'MIT', 'n_tokens_mistral': 5174, 'n_tokens_neox': 4886, 'n_words': 1530, 'path': 'Basic_Logics/Zillow/zillow.ipynb', 'repo_name': 'hshehjue/Machine_Learning_Projects', 'size': '15183'}