We have a data set of real estate prices from Aimes, Iowa. Our goal is to provide a model to predict prices using a subset of the columns to be able to reduce the size of the database the company is required to store. We are also interested in providing descriptive statistics for the market, specifically which neighborhoods have many sales and high prices. Finally, we are interested in describing how this data has changed over time.

Results summary

The most significant predictor of housing price is the OverallQuality of the property. Additionally, there is a nonlinear relation between this quality and SalePrice, with a quadratic term added to the model showing an even stronger correlation than the linear term.

I suggest a simplified model with six parameters to predict SalePrices. The linear model with ['LotArea', 'OverallQual', 'OverallCond', 'YearBuilt', 'GrLivArea', 'OverallQualSquared'] explains the variation in the data just about as well as the best model with twice as many parameters (difference in R-squared on a test set is less than 0.01). This parameter selection is stable over repeated selections of test and training data in a 2:1 split.

NAmes and CollgCr are the largest markets (with CollgCr also having high prices). NoRidge, NridgHt, and StoneBr have the highest average sale prices. Crawfor has prices that are trending higher.

Detailed Analysis

Read in the data

First save the list of variables we are allowed to use, and read in the just the training data from those columns.

avail_columns = ["LotArea", "Utilities", "Neighborhood", "BldgType", "HouseStyle", "OverallQual" \
           ,"OverallCond", "YearBuilt", "YearRemodAdd", "RoofStyle", "RoofMatl" \
           ,"GrLivArea", "FullBath", "HalfBath", "BedroomAbvGr", "KitchenAbvGr" \
           ,"MoSold", "YrSold", "SalePrice"]


df = pd.read_csv("train.csv", usecols=avail_columns)
print df

Explore each variable and correlations

Check the info on the data frame for column data types and null values. Any column with dtype of int64 is complete since a single string or NaN entry would change the dtype for the entire column. The object columns are all expected to have string values, so no worries there. Check the value_counts on the object columns for mis-spelled labels.

df['YrSold'].describe()

count    1460.000000
mean     2007.815753
std         1.328095
min      2006.000000
25%      2007.000000
50%      2008.000000
75%      2009.000000
max      2010.000000
Name: YrSold, dtype: float64

df.BldgType.value_counts()

1Fam      1220
TwnhsE     114
Duplex      52
Twnhs       43
2fmCon      31
Name: BldgType, dtype: int64

Need to decide which variables should become dummified. We do this with categorical classifications that we expect to have a strong impact on the result. Utilities and RoofMatl may not be good choices since there is so little variation… though we should inspect them to see if there is a strong correlation with the price.

correl = df.corr()**2
fig, ax = plt.subplots(figsize=(12,10))
sns.heatmap(correl);

for i in avail_columns:
    if df[i].dtype == int:
        df.plot(x=i,y="SalePrice",kind="scatter")
        plt.show()

Explore nonlinear relationships

Some of the features are correlated in a nonlinear fashion with the sales price, so I will check correlations with the logarithm of the sales price, and with the logarithm of the features… I can also try polynomial fits with the features. Overall Quality, YearBuilt,

Looks like some features show outliers (such as LotArea being greather than 60000 units, HalfBath ==2, or BedroomAbvGr > 5, KitchenAbvGr == 0 or 3). One can make a stron case that the estimate for these lots should be given by a different function than for the others… Perhaps these are times when the categorical variables would work well?

plt.scatter(df.OverallQual**2,df.SalePrice)
#df.corr()

<matplotlib.collections.PathCollection at 0x1178f5750>

df["OverallQualSquared"] = df.OverallQual**2
df["YearBuiltSquared"] = df.YearBuilt**2

Training test split

features = ['LotArea', 'OverallQual', 'OverallCond', 'YearBuilt', 'YearRemodAdd',
       'GrLivArea', 'FullBath', 'HalfBath',
       'BedroomAbvGr', 'KitchenAbvGr', 'MoSold', 'YrSold',
       'OverallQualSquared', 'YearBuiltSquared']

X = df[features]
y = df['SalePrice']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.33)

Fitting

# REGRESSION MODELS

results = []
for n in range(1,len(features)+1):
    max_r2 = 0
    best_features = []
    for cols in itertools.combinations(features, n):
        fit_xtrain = X_train[list(cols)]
        fit_xtest = X_test[list(cols)]        

#FIT LINEAR MODEL
        line = lm.LinearRegression(normalize=True)

#for algorithm in [line, lasso, ridge]:
        for algorithm in [line]:
            ## Fit model to data
            model = algorithm.fit(fit_xtrain,y_train)

            test_r2 = algorithm.score(fit_xtest, y_test)
#             print "Test r2 score : {}".format(test_r2);
            if test_r2 > max_r2:
                max_r2 = test_r2
                coeff_sum = abs(algorithm.coef_).mean()
                mse = mean_squared_error(y_true=y_test, y_pred=algorithm.predict(fit_xtest))
                best_features = list(cols)
    results.append( (max_r2,coeff_sum, mse, best_features))
results

[(0.64493256596024184,  3738.6815185367077,  2186715449.1880145,
  ['OverallQualSquared']),
 (0.72915774869635741,  1444.2671637627586,  1668006914.8561254,
  ['GrLivArea', 'OverallQualSquared']),
 (0.75483471919634171,   939.9441940699586,  1509872930.4412842,
  ['YearBuilt', 'GrLivArea', 'OverallQualSquared']),
 (0.79589913657403477,  12011.917742913709,  1256973939.1988173,
  ['LotArea',       'OverallQual',       'OverallCond',
   'YearBuilt',       'GrLivArea',       'HalfBath',
   'KitchenAbvGr',       'OverallQualSquared']),
 (0.79640356966477699,  10782.082461790924,  1253867341.6152034,
  ['LotArea',       'OverallQual',       'OverallCond',
   'YearBuilt',       'GrLivArea',       'HalfBath',
   'BedroomAbvGr',       'KitchenAbvGr',       'OverallQualSquared']),
 (0.79670268601380689,  9925.112054464651,  1252025206.0690444,
  ['LotArea',       'OverallQual',       'OverallCond',
   'YearBuilt',       'GrLivArea',       'FullBath',
   'HalfBath',       'BedroomAbvGr',       'KitchenAbvGr',
   'OverallQualSquared']),
 (0.79624488667288296,  7607.3176169180288,  1254844605.4153514,
  ['LotArea',       'OverallQual',       'OverallCond',
   'YearBuilt',       'YearRemodAdd',       'GrLivArea',
   'FullBath',       'HalfBath',       'BedroomAbvGr',
   'KitchenAbvGr',       'MoSold',       'YrSold',
   'OverallQualSquared',       'YearBuiltSquared'])]

Plot the progression of the best coefficient of determination as the number of features selected grows.

plt.plot(range(1,len(results)+1),[r[0] for r in results])
plt.xlabel("number of features")
plt.ylabel("best R^2")

Plot the absolute sum of the coefficients. Note that this only works when the features are normalized in some fashion before fitting. A significant increase in this sum is an indication of the model being overfit. Comparing this to the difference in R2 score would make a good cross check for overfitting.

plt.plot(range(1,len(results)+1),[r[1] for r in results])
plt.xlabel("number of features")
plt.ylabel("Abs sum of coeffs")

Plot the best mean squared error (MSE) for each number of features corresponding to the best R2 score.

plt.plot(range(1,len(results)+1),[r[2] for r in results])
plt.xlabel("number of features")
plt.ylabel("MSE")

Compare the be R2 score with a value from near the point of diminishing R2 returns.

print max([r[0] for r in results])
results[5]

0.796794359218


(0.79010957401318027,     12257.811452486383,     1292629492.7138808,
 ['LotArea',      'OverallQual',      'OverallCond',
  'YearBuilt',      'GrLivArea',      'OverallQualSquared'])

Selected model

The r-squared value for the test set drops plateaus around 0.83, when six features are used. The set selected is: [‘LotArea’, ‘OverallQual’, ‘OverallCond’, ‘YearBuilt’, ‘GrLivArea’, ‘OverallQualSquared’]. The linear model with these features selected from the entire training set will capture most of the trends in the data without risking overfitting or storing extraneous data (particualrly since one of the features selected is derived from another!).

The model with 13 parameters (excluding only month sold) consistently beats the model with 6 parameters, but only by a small margin. Considering the colinearity between many of the parameters (seen in the heatmap during initial data exploration) a smaller parameter set is desireable. Selecting parameters with Ridge and Lasso was not as successful as the fitting with limited parameter sets, so I have opted for the latter.

A model with 8 parameters sometimes did slightly better than the model with 6, but the feature selection was not consistent (due to colinearity), so I went for the six final parameters. The differences in R^2 between the full model and 6 parameters were someitmes as large as 0.01, while the 8 parameter fit was half that. This increase in success did not seem to justify the loss in certainty over the best feature set.

correl = df[['LotArea', 'OverallQual', 'OverallCond', 'YearBuilt', 'GrLivArea', 'OverallQualSquared','SalePrice']].corr()**2
fig, ax = plt.subplots(figsize=(12,10))
sns.heatmap(correl);

Neighborhood summary

The first figure shows the distribution of sales prices. Anything over $170k will be toward the high end of the market. Anything over $300k is really an exclusive property.

Below that is a box plot graph showing the distribution of prices in each neighborhood. It seems NoRidge, NridgHt, and StoneBr have higher median SalePrice than the other neighborhoods. A t test on this hypothesis shows significance with a p-value less than 0.0001.

Below that are tables showing the number of homes sold and the mean sale price in each neighborhood. Interestingly, there is no correlation between those two variables (as seen by the scatter plot below the tables.

NridgHt, in particular, has a good combination of sale price and sale volume. NAmes and CollgCr are two other neighborhoods that have lower home prices, but much more sales volume. CollgCr is a good market with over 100 homes sold in the dataset, and a mean sale price of nearly $200k.

Finally, I have some charts and tables showing how sale price and volume change over time in each neighborhood in the data set. This is a fairly stable market in terms of meaningful rankings of the neighborhoods. Crawfor, however, is worth watching, as it’s mean sale price has been trending upward during the duration of the data set.

To sum up: NAmes and CollgCr are the largest markets (with CollgCr also having high prices). NoRidge, NridgHt, and StoneBr have the highest average sale prices. Crawfor has prices that are trending higher.

df.SalePrice.hist(bins=20)

<matplotlib.axes._subplots.AxesSubplot at 0x1178356d0>

df.boxplot(column="SalePrice", by="Neighborhood",figsize=(14,8),rot=90,fontsize=16)

<matplotlib.axes._subplots.AxesSubplot at 0x1182aa550>

mask = (df['Neighborhood'] =='NoRidge' ) \
        | (df['Neighborhood']=='NridgHt') \
        | (df['Neighborhood']=='StoneBr')
tstat, pval = stats.ttest_ind(a= df[mask]['SalePrice'], \
                b= df[[not x for x in mask]]['SalePrice'], \
                equal_var=False)
print 2* pval

6.94365699791e-37

df.pivot_table(index=["Neighborhood"], values=["SalePrice"], \
               aggfunc=[np.mean,len]).sort_values(('mean','SalePrice'),ascending=False)

df.pivot_table(index=["Neighborhood"], values=["SalePrice"], \
               aggfunc=[np.mean,len]).sort_values(('len','SalePrice'),ascending=False)

df.pivot_table(index=["Neighborhood"], values=["SalePrice"], \
               aggfunc=[np.mean,len]).plot.scatter(x="mean",y="len")

<matplotlib.axes._subplots.AxesSubplot at 0x118173610>

df.pivot_table(index=["YrSold"], values=["SalePrice"], \
               columns=["Neighborhood"],aggfunc=np.mean \
              ).plot(figsize=(12,10))

<matplotlib.axes._subplots.AxesSubplot at 0x117bace50>

df.pivot_table(index=["Neighborhood"], values=["SalePrice"], \
               columns=["YrSold"],aggfunc=np.mean \
              ).sort_values(("SalePrice",2010),ascending=False)

df.pivot_table(index=["YrSold"], values=["SalePrice"], columns=["Neighborhood"],aggfunc=len \
              ).plot(figsize=(12,10))

<matplotlib.axes._subplots.AxesSubplot at 0x117da6f90>