<jupyter_start><jupyter_text>Project: Airbnb Seattle Analysis and PredictionIn this project, I wanggle and analyze the open data from Airbnb Seattle to to answer the following 4 questions: 1) How does the occupancy rate change during the record window?2) How does the price change during the record window?3) Are the trends of average price and occupancy rate the same?4) Using machine learning technique to predict the listing prices.<jupyter_code># import libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler 
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import LinearRegression
import seaborn as sns
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import r2_score, mean_squared_error
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Step 0: Load the DataThere are four files associated with this project (not including this one):- `listings.csv`- `calendar.csv`<jupyter_code>listings = pd.read_csv('listings.csv')
calendar = pd.read_csv('calendar.csv')<jupyter_output><empty_output><jupyter_text>Step 1: How does the occupancy rate change during the record window? Step 1.1 : PreprocessingThe infomation of the dataset **calendar** is checked. The dataframe  that describes the occupation state everyday is created from **calendar**. The data types of the features are corrected.<jupyter_code># Check the structure of the data 
calendar.info()
calendar.head()
calendar.iloc[-1,:]<jupyter_output><empty_output><jupyter_text>**Observation:** There are some missing values in the column 'price'. When the listing is not available (occupied), the price is NaN. The record window of the dataset is from 2016-01-04 to 2017-01-02. The avaiability of each listing is represented by strings 't' and 'f' instread of boolean variables. **Action:** Replace the strings by boolean variables. Since we care about the occupancy rate, the state should be converted from 'available' to 'occupied'.<jupyter_code># create new dataframe
occupy =  calendar[['date','available']]
occupy['occupied'] = occupy['available'].apply(lambda x: x== 'f')
occupy.drop('available',axis=1, inplace=True)

#change the data type of date
occupy.set_index('date')
occupy.index = pd.to_datetime(occupy.date)
occupy.drop('date',axis=1, inplace=True)
# check the dataframe
occupy.head()<jupyter_output><empty_output><jupyter_text>Step 1.2: The amount of listings in each dayBefore the analsis of occupancy rate, the number of listings in each day is checked.<jupyter_code>supply = occupy.resample('D').count()
plt.plot(supply['occupied']);
supply.head()<jupyter_output><empty_output><jupyter_text>**Observation:** The everyday number of listings is a constant. It implies that the dataset only contains a part of the data in seattle, which does not show the supply in the market.   Step 1.3: The trend of the occupancy rate<jupyter_code>#calculate the average occupancy rate everyweek
occupied_r = occupy.resample('W').mean()
#sketch the trend
ax1 = plt.plot(occupied_r['occupied'], color = 'b')
plt.ylabel('Occupancy rate', fontsize=12);<jupyter_output><empty_output><jupyter_text>**Observation:** The occupancy rate generally declines during 2016. It increases at the second quarter to a small extent but keeps dropping to the end of the year. The occupancy rate decreases from over 0.5 to around 0.25.   Step 2: How does the price change during the record window? Step 2.1 : PreprocessingAs mentioned in the last section, there are some missing values in price when the listing are not available. It makes sense that the price is not necessary to be showed when the listing is not available. Meanwhile, other available listings' prices indicate the market status. These missing values are all removed while the information is not lost. The original data type of price is not numeric which is converted in this section.<jupyter_code>#remove missing values
calendar_clean = calendar.dropna()
calendar_clean.drop(['available'],axis=1,inplace=True)

#convert data type of price to numeric
calendar_clean['price']=calendar_clean['price'].apply(lambda x: float(x.strip('$').replace(',','')))
calendar_clean['listing_id']=calendar_clean['listing_id'].astype(str)

#change the data type of date
calendar_clean.index = pd.to_datetime(calendar_clean.date)
calendar_clean.drop('date',axis = 1, inplace=True)

#check the conversion
calendar_clean.info()<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  This is separate from the ipykernel package so we can avoid doing imports until
C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:6: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  
C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:7: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/st[...]<jupyter_text>Step 2.2: The trend of the priceThe average price everyweek is calculated and sketched.<jupyter_code>#calculate the average price everyweek
price_mean = calendar_clean.resample('W').mean()

#sketch the trend
ax2=plt.plot(price_mean, color = 'r')
plt.ylabel('Price', fontsize=12);<jupyter_output><empty_output><jupyter_text>**Observation:** The average price first increases from around 120 to 150. Then, it fluctuates and falls suddenly after September. At the end of the year, the price is around 135.    Step 3: Are the trends of average price and occupancy rate the same?In this section, the plots in the above sections are merged for clear observation and comparison.<jupyter_code>fig, axes1= plt.subplots()
axes2 = axes1.twinx()
axes1.plot(occupied_r['occupied'],'b')
axes2.plot(price_mean, 'r')

axes1.set_ylabel('Ocuppancy rate', color='b');
axes2.set_ylabel('Price', color='r');

axes1.tick_params(axis='y', colors='b')
axes2.tick_params(axis='y', colors='r')<jupyter_output><empty_output><jupyter_text>**Observation:** Usually, when the occupancy rate is high, the renting cost is expensive. However, the trends of the occupancy rate and the price in Airbnb Seattle are not completely the same during 2016. In the first quater, the occupancy rate drops largely from over 0.5 to below 0.3, while the price increases with slight fluctuation. Then, The price keeps rising and reaches the peak in July. During this time, the occupancy rate slightly increases and remains stable from May to June then suddenly rises in July. For price, it remains the highest level around two months but sharply declines after September. The tendency of occupancy rate follows the similar step but it drops earlier. After the decrease, the price fluctuates at the level of 135. Although the price falls down, the occupancy rate does not recover. At the end of 2016, the occupancy rate drops to the bottom at 0.25.   Step 4: Predcting the price   Step 4.1: Preprocessing     step 4.1.1 remove missing values<jupyter_code>listings.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 3818 entries, 0 to 3817
Data columns (total 92 columns):
id                                  3818 non-null int64
listing_url                         3818 non-null object
scrape_id                           3818 non-null int64
last_scraped                        3818 non-null object
name                                3818 non-null object
summary                             3641 non-null object
space                               3249 non-null object
description                         3818 non-null object
experiences_offered                 3818 non-null object
neighborhood_overview               2786 non-null object
notes                               2212 non-null object
transit                             2884 non-null object
thumbnail_url                       3498 non-null object
medium_url                          3498 non-null object
picture_url                         3818 non-null object
xl_picture_url                      3498[...]<jupyter_text>**Obeservation:** Besides price, 91 features are recorded in this dataset. But not all of them are related to price.**Action:** Select the relevant features and create a new dataframe.<jupyter_code>#create new dataframe fro price predction
df = listings[['host_response_time', 'host_response_rate', 'host_acceptance_rate', 'host_is_superhost', 'neighbourhood_group_cleansed', 
               'property_type', 'room_type', 'accommodates','bathrooms', 'bedrooms','beds','bed_type', 'square_feet','price', 'guests_included', 'number_of_reviews',
              'review_scores_rating', 'review_scores_accuracy', 'review_scores_cleanliness', 'review_scores_checkin', 'review_scores_communication',
              'review_scores_location', 'review_scores_value', 'instant_bookable', 'cancellation_policy', 'reviews_per_month']]
def col_nan(df):
    """visualize the missing values propotion in each column.
    df -- dataframe
    """
    nan_prop = df.isnull().sum()/df.shape[0]
    plt.figure(figsize=(10,6));
    plt.bar(nan_prop.index,nan_prop);
    plt.xticks(rotation='vertical');
#check the proportion of the missing values in each column
col_nan(df)
def remove_nan(df):
    """remove the columns whose numbers of missing values are outliers.
    df -- dataframe 
    """
    nan_prop = df.isnull().sum()/df.shape[0]
    IQR = nan_prop.quantile(0.75) - nan_prop.quantile(0.25)
    #remove the columns of which the misiing value proportion is outliers.
    nan_prop = nan_prop[~(nan_prop > (nan_prop.quantile(0.75) + 1.5*IQR))]
    cols = list(nan_prop.index)
    df = df[df.columns.intersection(cols)]
    #remove the residual missing values
    df_wona = df.dropna()
    return df_wona
df_wona = remove_nan(df)
#check the processed dataframe
df_wona.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 2655 entries, 0 to 3814
Data columns (total 25 columns):
host_response_time              2655 non-null object
host_response_rate              2655 non-null object
host_acceptance_rate            2655 non-null object
host_is_superhost               2655 non-null object
neighbourhood_group_cleansed    2655 non-null object
property_type                   2655 non-null object
room_type                       2655 non-null object
accommodates                    2655 non-null int64
bathrooms                       2655 non-null float64
bedrooms                        2655 non-null float64
beds                            2655 non-null float64
bed_type                        2655 non-null object
price                           2655 non-null object
guests_included                 2655 non-null int64
number_of_reviews               2655 non-null int64
review_scores_rating            2655 non-null float64
review_scores_accuracy          2655 non-nul[...]<jupyter_text>**Observation:** All missing values have been removed. There are some nonnumerical features that should be further processed. However, some of them represent percentage, such as response rate and acceptance rate.  **Action:** Normalize the numerical features. Check the values of the nonnumerical features.   step 4.1.2 Feature normalization<jupyter_code>#select the numerical columns
num_cols = df_wona.select_dtypes(exclude=['object']).columns
df_num = df_wona[num_cols]
#visualize the distribution of these features
df_num.hist(figsize=(15,10));<jupyter_output><empty_output><jupyter_text>**Observation:** All the distributions are not Guassian distribution. Besides 'number_of_reviews', each feature contains the values of the same order of magnitude. **Action:** 'number_of_reviews' should be applied log-transfromation. Then, apply MinMaxScaler for normalization to all the numerical features.<jupyter_code>#log-transformation
df_num['number_of_reviews']=df_num['number_of_reviews'].apply(lambda x: np.log(x + 1))
df_num_scal = pd.DataFrame(data = df_num)
scaler = MinMaxScaler()
df_num_scal[df_num.columns] = scaler.fit_transform(df_num[df_num.columns])
df_num_scal.head()

df_wona.drop(num_cols,axis=1,inplace=True)
df_wona = pd.merge(df_wona,df_num_scal, left_index=True, right_index=True)<jupyter_output><empty_output><jupyter_text>step 4.1.3 Data types conversion<jupyter_code>#check the non-numerical values
non_num_cols = df_wona.select_dtypes(include=['object']).columns

for col in non_num_cols:
    print('{}: \n{}'.format(col,df_wona[col].unique()))<jupyter_output>host_response_time: 
['within a few hours' 'within an hour' 'within a day' 'a few days or more']
host_response_rate: 
['96%' '98%' '67%' '100%' '71%' '97%' '60%' '31%' '50%' '70%' '90%' '88%'
 '80%' '63%' '33%' '99%' '83%' '94%' '58%' '43%' '93%' '92%' '40%' '57%'
 '75%' '89%' '95%' '81%' '91%' '78%' '38%' '86%' '30%' '56%' '76%' '64%'
 '82%' '17%' '87%' '25%' '69%' '53%' '65%']
host_acceptance_rate: 
['100%' '0%']
host_is_superhost: 
['f' 't']
neighbourhood_group_cleansed: 
['<NAME>' 'Ballard' 'Other neighborhoods' 'Cascade' 'Central Area'
 'University District' 'Downtown' 'Magnolia' 'West Seattle' 'Interbay'
 'Beacon Hill' 'Rainier Valley' 'Delridge' 'Seward Park' 'Northgate'
 'Capitol Hill' 'Lake City']
property_type: 
['Apartment' 'House' 'Cabin' 'Condominium' 'Camper/RV' 'Bungalow'
 'Townhouse' 'Loft' 'Boat' 'Other' 'Dorm' 'Bed & Breakfast' 'Yurt'
 'Chalet' 'Tent' 'Treehouse']
room_type: 
['Entire home/apt' 'Private room' 'Shared room']
bed_type: 
['Real Bed' 'Futon' 'Pull-out[...]<jupyter_text>**Observation:** Some features represent percentage, such as response rate and acceptance rate, which should be converted to numerical values. Some features are represented by strings 't' and 'f' instread of boolean variables. **Action:** Converte the feature types.<jupyter_code>def non_num_conv(non_num_cols, df_wona):
    """
    non_num_cols: nonnumerical features names (list)
    df_wona: dataframe to be processed
    """
    
    rate_cols = [col for col in non_num_cols if 'rate' in col]
    for rate_col in rate_cols:
        df_wona[rate_col]=df_wona[rate_col].apply(lambda x: float(x.strip('%'))/100)
    
    tf_cols = [col for col in non_num_cols if 't' in df_wona[col].unique()]

    for tf_col in tf_cols:
    
        df_wona[tf_col]=df_wona[tf_col].apply(lambda x: x == 't')
    
    df_wona['price'] = df_wona['price'].apply(lambda x: float(x.strip('$')))
    return df_wona
df_wona = non_num_conv(non_num_cols, df_wona)

#check converted results
for col in non_num_cols:
    print('{}: \n{}'.format(col,df_wona[col].unique()))<jupyter_output>host_response_time: 
['within a few hours' 'within an hour' 'within a day' 'a few days or more']
host_response_rate: 
[0.96 0.98 0.67 1.   0.71 0.97 0.6  0.31 0.5  0.7  0.9  0.88 0.8  0.63
 0.33 0.99 0.83 0.94 0.58 0.43 0.93 0.92 0.4  0.57 0.75 0.89 0.95 0.81
 0.91 0.78 0.38 0.86 0.3  0.56 0.76 0.64 0.82 0.17 0.87 0.25 0.69 0.53
 0.65]
host_acceptance_rate: 
[1. 0.]
host_is_superhost: 
[False  True]
neighbourhood_group_cleansed: 
['<NAME>' 'Ballard' 'Other neighborhoods' 'Cascade' 'Central Area'
 'University District' 'Downtown' 'Magnolia' 'West Seattle' 'Interbay'
 'Beacon Hill' 'Rainier Valley' 'Delridge' 'Seward Park' 'Northgate'
 'Capitol Hill' 'Lake City']
property_type: 
['Apartment' 'House' 'Cabin' 'Condominium' 'Camper/RV' 'Bungalow'
 'Townhouse' 'Loft' 'Boat' 'Other' 'Dorm' 'Bed & Breakfast' 'Yurt'
 'Chalet' 'Tent' 'Treehouse']
room_type: 
['Entire home/apt' 'Private room' 'Shared room']
bed_type: 
['Real Bed' 'Futon' 'Pull-out Sofa' 'Airbed' 'Couch']
price: 
[ 85. 150. 97[...]<jupyter_text>**Observation:** There are some categorical features**Action:** Apply one-hot ecoding<jupyter_code>#update non_num columns list
non_num_cols = df_wona.select_dtypes(include=['object']).columns
#get dummy variables
df_non_num = df_wona[non_num_cols]
dummy_vars = pd.get_dummies(df_non_num)
df_wona.drop(non_num_cols,axis=1,inplace=True)
df_cleaned=pd.merge(df_wona,dummy_vars, left_index=True, right_index=True)

# check the converted results
df_cleaned.head()
df_cleaned.describe()<jupyter_output><empty_output><jupyter_text>Step 4.2: Price distribution visualization<jupyter_code>df_cleaned['price'].hist(bins=50);
plt.ylabel('Count');
plt.xlabel('Listing price in $');
plt.title('Histogram of listing prices');<jupyter_output><empty_output><jupyter_text>**Observation:** It is a right skewed distribution. Most prices are within the range from 0 to 200. The most expensive listings almost reach 1000.<jupyter_code># mode of price
df_cleaned['price'].mode()
q_outlier = 0.99
df_cleaned['price'].quantile(q_outlier)
df_4_reg = df_cleaned[df_cleaned['price']<df_cleaned['price'].quantile(q_outlier)]
df_4_reg['price'].hist(bins=50);<jupyter_output><empty_output><jupyter_text>Step 4.3: Regression by machine learningIn this section, two kinds of regression algorithm are applied, linear regression model and random forest model. The results are evaluated by R2 score and mean square error.<jupyter_code>X = df_4_reg.drop('price', axis=1, inplace=False)
y = df_4_reg['price'].reshape(-1,1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:2: FutureWarning: reshape is deprecated and will raise in a subsequent release. Please use .values.reshape(...) instead<jupyter_text>Step 4.3.1 Linear regression model<jupyter_code>lm_model = LinearRegression()
lm_model.fit(X_train,y_train)
y_train_pred = lm_model.predict(X_train)
y_test_pred = lm_model.predict(X_test)
print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))<jupyter_output>MSE train: 1315.307, test: 55733161759341867224268800.000
R^2 train: 0.688, test: -12414568999121368121344.000<jupyter_text>**Observation:** the linear regression model did not work well.   Step 4.3.2 Random forest model<jupyter_code>Rforest = RandomForestRegressor(n_estimators=200, 
                               random_state=1,
                               n_jobs=-1)
Rforest.fit(X_train, y_train)
y_train_pred = Rforest.predict(X_train).reshape(-1,1)
y_test_pred = Rforest.predict(X_test).reshape(-1,1)
print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))<jupyter_output>MSE train: 170.646, test: 1546.774
R^2 train: 0.959, test: 0.655<jupyter_text>Step 4.3.3 Parameters tuningThe random forest model is selected for parameters tuning.<jupyter_code>param_grid={'n_estimators':[50,100,200],
            'max_depth':[25,30,35], 
            'min_samples_leaf':[3,5], 
            #'max_features':[0.5],
            'random_state':[1],
            'n_jobs':[-1]}
grid_obj = GridSearchCV(RandomForestRegressor(),param_grid, scoring = 'r2')
grid_obj.fit(X_train, y_train)
best_r = grid_obj.best_estimator_
best_predictions = best_r.predict(X_test).reshape(-1,1)
y_train_grid = best_r.predict(X_train).reshape(-1,1)
print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_grid),
        mean_squared_error(y_test, best_predictions)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_grid),
        r2_score(y_test, best_predictions)))
    
print (grid_obj.best_estimator_)
print (Rforest)
plt.scatter(y_train_grid, y_train_grid - y_train,
            c='blue', marker='o', label='Training data')
plt.scatter(best_predictions,  best_predictions - y_test,
            c='lightgreen', marker='s', label='Test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')<jupyter_output><empty_output><jupyter_text>Step 4.3.4 Feature importancesIn this part, the two-5 importances of features when making predictions are ranked.<jupyter_code># feature importances
importances = best_r.feature_importances_
def feature_plot(importances, X_train):
    """plot the feature importances of the top-5 most weighted features
    importances -- output of feature importances 
    X_train -- training dataset without label
    """
    
    # Display the five most important features
    indices = np.argsort(importances)[::-1]
    columns = X_train.columns.values[indices[:5]]
    values = importances[indices][:5]

    # Creat the plot
    fig = plt.figure(figsize = (9,6))
    plt.title("Normalized Weights for First Five Most Predictive Features", fontsize = 16)
    plt.bar(np.arange(5), values, width = 0.6, align="center", color = '#00A000', \
          label = "Feature Weight")
    plt.bar(np.arange(5) - 0.3, np.cumsum(values), width = 0.2, align = "center", color = '#00A0A0', \
          label = "Cumulative Feature Weight")
    plt.xticks(np.arange(5), columns, rotation=20)
    plt.xlim((-0.5, 4.5))
    plt.ylabel("Weight", fontsize = 12)
    plt.xlabel("Feature", fontsize = 12)
    
    plt.legend(loc = 'upper left')
    plt.tight_layout()
    plt.show()  
feature_plot(importances, X_train)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11218524144915527, 'fraction_numerical': 0.03489389126732035, 'mean_word_length': 4.447342026078235, 'pattern_counts': {'":': 0, '<': 90, '<?xml version=': 0, '>': 90, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '597203', 'license': 'CNRI-Python', 'n_tokens_mistral': 7638, 'n_tokens_neox': 7014, 'n_words': 2039, 'path': 'Seattle Airbnb.ipynb', 'repo_name': 'JaneRu1123/Airbnb-Seattle', 'size': '21731'}