<jupyter_start><jupyter_text>Toy Tweet Project_<NAME>, August 2018_<jupyter_code>import time
notebookstart= time.time()

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import os
import gc
from IPython.display import display
print("Data:\n",os.listdir("../input/tweets/"))

# Models Packages
from sklearn import metrics
from sklearn.metrics import mean_squared_error
from sklearn import feature_selection
from sklearn.model_selection import train_test_split
from sklearn import preprocessing
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error, make_scorer
from sklearn.ensemble import RandomForestRegressor
# Gradient Boosting
import lightgbm as lgb

# Tf-Idf
from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer
from sklearn.pipeline import FeatureUnion
from scipy.sparse import hstack, csr_matrix
from nltk.corpus import stopwords 

# Visualization
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline

# Ignore Future Warnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
# Load Data
PATH = "../input/tweets/" # Google's Cloud Directory (Linux Backend)
df = pd.read_excel(PATH + "tweets.xlsx", sheet_name="Sheet1")#.sample(5000) # Uncomment to Debug..
print("Dataset Shape: {} Rows, {} Columns".format(*df.shape))

# Parse Date Time
df["TweetPostedTime"] = pd.to_datetime(df['TweetPostedTime'])
df["UserSignupDate"] = pd.to_datetime(df['UserSignupDate'])

# Glance
display(df.head())<jupyter_output><empty_output><jupyter_text>Pre-Processing..1. Missing Values2. Remove Retweets<jupyter_code># Remove Columns with 95%+ Missing
missing = round(df.isnull().sum()/ df.shape[0]*100).reset_index().rename({"index":"columns",0:"missing"}, axis =1 )
high_missing_columns = missing.loc[missing.missing > 95, "columns"]
print("Columns to remove (95% missing Values and Over)\n", list(high_missing_columns))
df.drop(high_missing_columns,axis =1, inplace= True)
# Remove Retweets..
print("Retweet Percent Breakdown:\n", round(df["TweetRetweetFlag"].value_counts(normalize=True)*100))
df = df.loc[df["TweetRetweetFlag"] == False, :]
print("\nDataframe Dimension after removing retweets: {} Rows, {} Columns".format(*df.shape))
df.drop("TweetRetweetFlag", axis = 1, inplace=True) # Obsolete variable
def custom_describe(df):
    """
    Custom Describe function to extract important information
    """
    unique_count = []
    for x in df.columns:
        mode = df[x].mode().iloc[0]
        unique_count.append([x,
                             len(df[x].unique()),
                             df[x].isnull().sum(),
                             mode,
                             df[x][df[x]==mode].count(),
                             df[x].dtypes])
    print("Dataframe Dimension: {} Rows, {} Columns".format(*df.shape))
    return pd.DataFrame(unique_count, columns=["Column","Number of Unique Values","Missing Count","Mode","Mode Occurence","dtype"]).set_index("Column")

# Visualize
custom_describe(df)<jupyter_output><empty_output><jupyter_text>Feature Engineering:1. Time Features and Account Age2. Extract URL Domain3. Label Encoding and Removing Redundant Variables.. 4. Count Vectorizing each Hastag in the TweetHastag feature5. Meta Text Features<jupyter_code># Source: My Taxi Far Prediction Project - https://www.kaggle.com/nicapotato/taxi-rides-time-analysis-and-oof-lgbm
def prepare_time_features(df, time_var):
    df['hour_of_day'] = df[time_var].dt.hour
    df['week'] = df[time_var].dt.week
    df['month'] = df[time_var].dt.month
    df["year"] = df[time_var].dt.year
    df['day_of_year'] = df[time_var].dt.dayofyear
    df['week_of_year'] = df[time_var].dt.weekofyear
    df["weekday"] = df[time_var].dt.weekday
    df["quarter"] = df[time_var].dt.quarter
    df["day_of_month"] = df[time_var].dt.day
    
    return df
# Run Function..
df = prepare_time_features(df, "TweetPostedTime")

# Calculate Account Age- The age is relative to most recent tweet in df
df["account_age_days"] = (df.TweetPostedTime.max()- df["UserSignupDate"]).astype('timedelta64[D]').astype(int)
# Remove time-formatted variables
df.drop(['TweetPostedTime', "UserSignupDate"], axis = 1, inplace=True)
# Extract Tweet Source
df["TweetSource"] = df["TweetSource"].apply(lambda url: url.split(">")[1].split("</a")[0])
print(df["TweetSource"].value_counts()[:5])
# Drop Redundant Variables
df.drop(["UserName","UserScreenName", "TweetID", "UserLink"], axis =1, inplace=True)

# Label Encoder:
categorical_to_encode = ["UserID", "UserExpandedLink", "TweetSource"]
lbl = preprocessing.LabelEncoder()
for col in categorical_to_encode:
    df[col] = lbl.fit_transform(df[col].astype(str))
# Seperating Features for Modeling
text_features = ["UserLocation","TweetBody", "UserDescription"]
basic_features = [x for x in df.columns if x not in text_features + ["TweetHashtags", "TweetID"]]
# Efficient One Hot Encoder for Hastags
# Source: My Cuisine Predicting Project - https://www.kaggle.com/nicapotato/this-model-is-bland-simple-logistic-starter
vect = CountVectorizer(tokenizer=lambda x: [i.strip() for i in x.split(',')], lowercase=False)
dummies = vect.fit_transform(df['TweetHashtags'].str.lower().apply(lambda x: [str(x)]).apply(','.join))
print("Vectorized Hashtag Sample : ", vect.get_feature_names()[14:22])
vect_names = vect.get_feature_names()

# Merge new features and Delete original variable
df = pd.concat([df.reset_index(drop=True), pd.DataFrame(dummies.todense(),columns=vect.get_feature_names())], axis=1)
print("Vocab Length: ", len(vect.get_feature_names()))
print("All Data Shape: ", df.shape)
df.drop("TweetHashtags",axis=1,inplace=True)
del dummies; gc.collect();
# Once again my work - Avito NLP Model
# https://www.kaggle.com/nicapotato/bow-meta-text-and-dense-features-lgbm/code
for cols in text_features:
    df[cols] = df[cols].astype(str) 
    df[cols] = df[cols].astype(str).fillna('FILLNA') # FILL NA
    df[cols] = df[cols].str.lower() # Lowercase all text, so that capitalized words dont get treated differently
    df[cols + '_num_words'] = df[cols].apply(lambda comment: len(comment.split())) # Count number of Words
    df[cols + '_num_unique_words'] = df[cols].apply(lambda comment: len(set(w for w in comment.split())))
    df[cols + '_words_vs_unique'] = df[cols+'_num_unique_words'] / df[cols+'_num_words'] * 100 # Count Unique Words<jupyter_output><empty_output><jupyter_text>Some Visualizations1. Dependent Variable Distribution2. Correlations<jupyter_code>f, ax = plt.subplots(figsize = [8,4])
sns.kdeplot(df["TweetRetweetCount"], ax= ax)
ax.set_title("RetweetCount Distribution (Continuous Dependent Variable)")
ax.set_xlabel("Retweets")
ax.set_ylabel("Density Distribution")
plt.show()<jupyter_output><empty_output><jupyter_text>**Correlations:**<jupyter_code>f, ax = plt.subplots(figsize=[10,7])
date_feats = ["week","month","hour_of_day","year","day_of_year","week_of_year","weekday","quarter","day_of_month"]
sns.heatmap(df[[x for x in basic_features if x not in date_feats + ["UserID","TweetSource"]]].corr(),
            annot=False, fmt=".2f",cbar_kws={'label': 'Correlation Coefficient'},cmap="plasma",ax=ax, linewidths=.5)
ax.set_title("Dense Features Correlation Matrix")
plt.savefig('correlation_matrix.png')<jupyter_output><empty_output><jupyter_text>**Red Flag - TweetFavoriteCount may be Leaky**Reflection: Tweet Favorite Count is similar to Tweet Retweet Count in the sense that they represent the response of the twitter community after the tweet is posted.  My understanding of this problem is to predict the retweet potential of the tweet characteristics in itself. Therefore, I shall **remove the TweetFavoritesCount** variable. Note that if this were more of a inferential project, I would keep the variable and measure the effect of the causation and feedback loop between these variables.<jupyter_code>print("Correlation Coefficient: ",df["TweetRetweetCount"].corr(df["TweetFavoritesCount"]))
f, ax = plt.subplots(figsize = [8,4])
sns.regplot(data = df, x = "TweetRetweetCount", y = "TweetFavoritesCount", ax= ax)
ax.set_title("Retweet Count and Favorite Count Correlation..")
plt.show()
df.drop(["TweetFavoritesCount"], axis = 1, inplace = True)<jupyter_output><empty_output><jupyter_text>Modeling without Text1. Random Forest2. Gradient Boosting<jupyter_code># Mean Square Error - Model Selection Metric
def mse(y_true, y_pred):
    return np.mean((y_true-y_pred)**2)
scoring = make_scorer(mse, greater_is_better=False)

# Helper Function to visualize feature importance for my Decision Trees..
# Source: My Titanic ML guide - https://www.kaggle.com/nicapotato/titanic-voting-pipeline-stack-and-guide
plt.rcParams['figure.figsize'] = (8, 8)
def test_score_and_feature_imp(model, predictors):
    MO = model.fit(X_train, y_train)
    print("Test Score: ", mse(y_test, model.predict(X_test)))
    feat_imp = pd.DataFrame({"Features" : predictors, "Importance" : MO.feature_importances_}).sort_values(by = "Importance",ascending=False)[:30]
    sns.barplot(x="Importance", y="Features", data=feat_imp)
    plt.title("Feature Importance")
    plt.show()<jupyter_output><empty_output><jupyter_text>Simple Basline: Random Forest<jupyter_code>print("Missing Values?",df[[x for x in df.columns if x not in text_features + ["TweetHashtags", "TweetID"]]].isnull().sum().any())
features = [x for x in df.columns if x not in text_features + ["TweetHashtags", "TweetID"]]

# Define y/x
y = df["TweetRetweetCount"]
X = df[[x for x in features if x not in ["TweetRetweetCount"]]]

# Train / Test Split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.3, random_state=11)

print("Train/Test - X/y Shape")
X_train.shape, y_train.shape, X_test.shape, y_test.shape
regressor = RandomForestRegressor(random_state=0, n_estimators = 10, max_depth = 4)
score = cross_val_score(regressor, X, y, cv=5, scoring = scoring)
print("CV - Mean: {} +/- {}".format(round(-score.mean(), 2), round(score.std(),2 )))
test_score_and_feature_imp(regressor, X.columns)<jupyter_output><empty_output><jupyter_text>Randomforest seems to put all of the prediction power into UserFollersCount, furthermore, it experiences massive performance variance accross cross-validation folds. **Light Gradient Boosting Model**<jupyter_code># Source: My Seed Diversification Notebook for Santander Bank
# https://www.kaggle.com/nicapotato/lgbm-cv-tuning-and-seed-diversification
lgtrain = lgb.Dataset(df[features],y ,feature_name = "auto")

print("Light Gradient Boosting Regressor: ")
lgbm_params =  {
    'task': 'train',
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': 'mse',
                }

modelstart= time.time()
# Find Optimal Parameters / Boosting Rounds
lgb_cv = lgb.cv(
    params = lgbm_params,
    train_set = lgtrain,
    num_boost_round=2000,
    stratified=False,
    nfold = 5,
    verbose_eval=150,
    seed = 23,
    early_stopping_rounds=75)

mymetric = "l2"
optimal_rounds = np.argmin(lgb_cv[mymetric + '-mean'])
best_cv_score = min(lgb_cv[mymetric + '-mean'])

print("\nOptimal Round: {}\nOptimal Score: {} + {}".format(
    optimal_rounds,best_cv_score,lgb_cv[mymetric + '-stdv'][optimal_rounds]))
print("Model Runtime: %0.2f Minutes"%((time.time() - modelstart)/60))
lgb_clf = lgb.train(
    lgbm_params,
    lgb.Dataset(X_train,y_train ,feature_name = "auto"),
    num_boost_round=optimal_rounds + 1,
    verbose_eval=100
)
print("Test Score: ", mse(y_test, lgb_clf.predict(X_test)))

# Feature Importance Plot
f, ax = plt.subplots(figsize=[7,10])
lgb.plot_importance(lgb_clf, max_num_features=50, ax=ax)
plt.title("Light GBM Feature Importance")
plt.show()<jupyter_output><empty_output><jupyter_text>*** Text Processing:1. Process Tweet Body, User Description with TF-IDF and user location with CountVectorizer.    - Character grams chosen to capture different languages and twitter slang lingo<jupyter_code># Free Up Memory with Garbage Collector..
del X, X_train, X_test, y_train, y_test, regressor; gc.collect()
# Source:
# My Advertisement Quality Prediction Algorithm for Avito
# https://www.kaggle.com/nicapotato/bow-meta-text-and-dense-features-lgbm/code

# # Prep Text
text_df = df.loc[:,text_features]
    
# TF-IDF
default_preprocessor = CountVectorizer().build_preprocessor()
def build_preprocessor(field):
    field_idx = list(text_df.columns).index(field)
    return lambda x: default_preprocessor(x[field_idx])

tfidf_para = {
    "stop_words": "english",
    "analyzer": 'char',
    "token_pattern": r'\w{1,}',
    "ngram_range": (2, 6),
    "sublinear_tf": True,
    "dtype": np.float32,
    "norm": 'l2',
    "smooth_idf":False
}

vectorizer = FeatureUnion([
        ('TweetBody',TfidfVectorizer(
            max_features= 500,
            **tfidf_para,
            preprocessor=build_preprocessor('TweetBody'))),
        ('UserDescription',TfidfVectorizer(
            **tfidf_para,
            max_features=500,
            preprocessor=build_preprocessor('UserDescription'))),
        ('UserLocation',CountVectorizer(
            ngram_range=(1, 1),
            max_features=100,
            preprocessor=build_preprocessor('UserLocation')))
    ])

start_vect=time.time()
text_df = vectorizer.fit_transform(text_df.values)
text_vocab = vectorizer.get_feature_names()
print("Vectorization Runtime: %0.2f Minutes"%((time.time() - start_vect)/60))
# Remove Raw Text Features
df.drop(text_features + ["TweetRetweetCount"], axis=1 ,inplace= True)

# Sparse Matrix
fulldf = hstack([csr_matrix(df.values), text_df])
fulldf_feature_names = list(df.columns) + text_vocab
print("Final DF Shape: {} Rows, {} Columns".format(*fulldf.shape))
print("Structured Features: ", len(df.columns))
print("Text Features: ", len(text_vocab))

# Clear Memory
del df, text_df; gc.collect();<jupyter_output><empty_output><jupyter_text>All Features Model1. Gradient Boosting - Model chosen because of its ability to efficiently process large/wide data.<jupyter_code># Source: My Seed Diversification Notebook for Santander Bank
# https://www.kaggle.com/nicapotato/lgbm-cv-tuning-and-seed-diversification
lgtrain = lgb.Dataset(fulldf,y ,feature_name = fulldf_feature_names)

print("Light Gradient Boosting Regressor: ")
lgbm_params =  {
    'task': 'train',
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': 'mse',
                }

modelstart= time.time()
# Find Optimal Parameters / Boosting Rounds
lgb_cv = lgb.cv(
    params = lgbm_params,
    train_set = lgtrain,
    num_boost_round=2000,
    stratified=False,
    nfold=5,
    verbose_eval=150,
    seed=23,
    early_stopping_rounds=75)

mymetric = "l2"
optimal_rounds = np.argmin(lgb_cv[mymetric + '-mean'])
best_cv_score = min(lgb_cv[mymetric + '-mean'])

print("\nOptimal Round: {}\nOptimal Score: {} + {}".format(
    optimal_rounds,best_cv_score,lgb_cv[mymetric + '-stdv'][optimal_rounds]))
print("Model Runtime: %0.2f Minutes"%((time.time() - modelstart)/60))
# Test / Train / Split
X_train, X_test, y_train, y_test = train_test_split(fulldf, y, test_size=.3, random_state=11)

# Run on Train Set
lgb_clf = lgb.train(
    lgbm_params,
    lgb.Dataset(X_train,y_train ,feature_name = fulldf_feature_names),
    num_boost_round=optimal_rounds + 1,
    verbose_eval=100
)
# Test Score
print("Test Score: ", mse(y_test, lgb_clf.predict(X_test)))

# Feature Importance Plot
f, ax = plt.subplots(figsize=[7,10])
lgb.plot_importance(lgb_clf, max_num_features=50, ax=ax)
plt.title("Light GBM Feature Importance")
plt.show()
print("Notebook Runtime: %0.2f Minutes"%((time.time() - notebookstart)/60))<jupyter_output><empty_output><jupyter_text>***  Toy Analyst Task¶_<NAME>, August 2018_**Assignment: **1. Can you identify which service has the most influence on growth in each sector from 2015-2019?2. Which firm type is having the biggest influence on growth in the top 5 fastest growing countries between 2015-2019?3. Can you forecast the current data set out to 2022? Use any methodology and external data sources you think are relevant.4. Present these results in whatever way you think is best, eg; dashboard, charts, graphs, tables, etc , If you can show your proficiency in R, and Tableau, that would be great but feel free to use any other programme you are comfortable with.<jupyter_code># Math and Data Processing
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import random
import gc
import os
import time
notebookstart= time.time()

# Visualization
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
from IPython.display import display

# Ignore Future Warnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
print("Data:\n",os.listdir("../input/sgr-analyst"))
# Load Data into Memory
df = pd.read_excel("../input/sgr-analyst/SGR Analyst Task.xlsx", sheet_name="Consulting market size")
print("Dataset:")
display(df.head())<jupyter_output><empty_output><jupyter_text>Quick Preliminary Data ExplorationThe first thing I want to do is explore the distribution for my features. This will enable me to understand the consistency of the dataset.Going forward, I am assuming that the yearly values represent the yearly growth in millions of dollars.<jupyter_code>f, ax = plt.subplots(1,2,figsize=[12,5])
# Figure 1
sns.countplot(y = df["Region"], order= df["Region"].value_counts().index[:], ax=ax[0])
ax[0].set_title("Region Occurence Count")
# Figure 2
sns.countplot(y = df["Report region"], order= df["Report region"].value_counts().index[:], ax=ax[1])
ax[1].set_title("Report Region Occurence Count")
# Format
plt.tight_layout(pad=0)
plt.show()<jupyter_output><empty_output><jupyter_text>**Additional Information About Categorical Variables:** The rest of the categorical variables Country, Industry, Sector, Service, and Firm Type are equally distributed across their respective classes. Next I want to see what their subclasses are.<jupyter_code>print("Glancing at Classes..")
for x in ["Country","Industry","Sector","Service","Firm Type"]:
    print("Classes in the {} Variable:\n {}\n".format(x,df[x].unique()))<jupyter_output><empty_output><jupyter_text>**Glancing at the Growth Rates** But first, I must stack the data so that the growth variable has one column and the year variable has its own column.<jupyter_code># Rename Variables
df.columns = ['Region', 'Report region', 'Country', 'Industry', 'Sector', 'Service',
       'Firm Type', '2013', '2014', '2015', '2016', '2017',
       '2018', '2019']

# Stack Growth Variable into a single variable
melt_df = pd.melt(df, id_vars=list(df.columns[:7]), value_vars=list(df.columns[7:]), value_name='Growth', var_name = "Year")
melt_df["Year"] = melt_df["Year"].astype(int)
print("Stack Data")
display(melt_df.head())<jupyter_output><empty_output><jupyter_text>**Understanding the Growth Variable:** These lines are very close together. It appears like as the years increase, the kurtosis decreases and the right tail becomes more extreme. This suggests growth increases over time, exponentially.<jupyter_code># Density Plot by Year
f, ax = plt.subplots(figsize=[12,5])
for year in melt_df.Year.unique():
    sns.kdeplot(melt_df.loc[melt_df.Year == year, "Growth"], label = year)
ax.set_title("Growth Density Distribution")
ax.set_xlabel("Growth")
ax.set_ylabel("Density")
plt.show()<jupyter_output><empty_output><jupyter_text>**Growth and Categorical Variables:** I am curious to see how growth interacts with variables such as Region, Industry, Service, and Firm Type.<jupyter_code>plot_vars = ["Region","Industry","Service","Firm Type"]
f, ax = plt.subplots(len(plot_vars)+ 1,1, figsize = [10,16], sharey = True)

# All Data
sns.pointplot(x="Year", y="Growth",data=melt_df, err_style="bars", ci=68, ax=ax[0])
ax[0].set_title("Growth Average by Year")

# By Specific Variable
for i,var in enumerate(plot_vars):
    sns.pointplot(x= "Year", y= "Growth", hue= var, data= melt_df, err_style="bars", ci= 67, ax= ax[i+1])
    ax[i+1].set_title("{} Growth Average by Year".format(var.capitalize()))
    ax[i+1].legend(fontsize='large', loc='center left', bbox_to_anchor=(1, 0.5))
plt.tight_layout(pad=0)
plt.show()<jupyter_output><empty_output><jupyter_text>**Interpretation:** In general, the trends points to a slight growth acceleration, with classes of diffrent growth levels and growth acceleration patterns.*** Q1: Can you identify which service has the most influence on growth in each sector from 2015-2019?¶Judging by the Service Growth Average by Year plot above, it appears like Risk is the dominant service across the board, although HR and Technology have seen a rate of accelerated growth from 2017 and onwards, so the trend suggests that they may overtake the impact of Risk in the next couple years.**Create New Features: Growth Difference Year by Year** Snapshots of dataset progression to facilitate explanation.<jupyter_code>print("Original Dataset")
display(df.head())
diff = df.iloc[:,-7:].diff(axis=1).iloc[:,1:]

print("Measuring Difference Year to Year Dataset")
diff = pd.concat([df.iloc[:,:7], diff], axis=1)
display(diff.head())

print("Stack Growth Difference Dataset")
diff_melt_df = pd.melt(diff, id_vars=list(diff.columns[:7]), value_vars=list(diff.columns[7:]), value_name='Growth', var_name = "Year")
diff_melt_df["Year"] = diff_melt_df["Year"].astype(int)
display(diff_melt_df.head())<jupyter_output><empty_output><jupyter_text>**Glancing at Sector Growth: Defining Influence in Two Different Ways**  Since there are 29 Sectors, there will be too many lines to attach to a name. My first step is to visualize all the lines to get a overview.<jupyter_code>f, ax = plt.subplots(1,2, figsize = [14,5], sharex= False)
sns.pointplot(x="Year", y="Growth", hue= "Sector", data=melt_df, ax=ax[0], ci=None)
ax[0].set_title("Sector Growth Average by Year")
ax[0].legend_.remove()

sns.pointplot(x="Year", y="Growth", hue= "Sector", data = diff_melt_df, ax=ax[1], ci=None)
ax[1].set_title("Sector Average Growth Difference by Year")
ax[1].set_ylabel("Growth Difference")
ax[1].legend_.remove()

plt.show()<jupyter_output><empty_output><jupyter_text>**Interpretation: ** The plot above serves to understand the general trend, as well as to spot potential outliers. An interesting finding here is that growth rates differ between sectors. Furthermore, there appears to be some sectors that experience low growth for a few years before sky rocketing. This is an important consideration when ranking sectors because some metrics may not effectively capture this.**Action:** Therefore, I will create a new feature which tracks growth change over time periods in order to the capture up and coming aspects of the data [Plot 2]. To capture general influence over the years, I will use the arithemtic mean [Plot 1].**Result: ** As I suspected, when ranking by average growth difference, different sectors appear at the very top, and there is a even more drastic difference between sectors.**Heatmap to Weigh Service Influence: **Heatmaps are my go-to visualization for variables with 10+ categories.<jupyter_code># Compute Average Growth
df["mean"] = df.loc[:,['2013', '2014', '2015', '2016', '2017', '2018', '2019']].mean(axis = 1)
df["stdev"] =  df.loc[:,['2013', '2014', '2015', '2016', '2017', '2018', '2019']].std(axis = 1)

# Compute Average Growth Difference
diff["mean"] =  diff.loc[:,['2014', '2015', '2016', '2017', '2018', '2019']].mean(axis = 1)
diff["stdev"] =  diff.loc[:,['2014', '2015', '2016', '2017', '2018', '2019']].std(axis = 1)

# Heatmap
f, ax = plt.subplots(1,2, figsize=[15,7], sharey=True)
# Plot 1
tps = df.pivot_table(values=['mean'], index='Sector',columns='Service',aggfunc='mean')
tps.columns = tps.columns.droplevel()
tps = tps.div(tps.sum(axis=1), axis=0)
sns.heatmap(tps, annot=False, fmt=".2f",cmap= sns.light_palette("blue", as_cmap=True), ax=ax[0])
ax[0].set_title("Average Growth's\nInfluence of Services on Sectors")
del tps
# Plot 2
tps = diff.pivot_table(values=['mean'], index='Sector',columns='Service',aggfunc='mean')
tps.columns = tps.columns.droplevel()
tps = tps.div(tps.sum(axis=1), axis=0)
sns.heatmap(tps, annot=False, fmt=".2f",cbar_kws={'label': 'Percent Influence'},cmap= sns.light_palette("blue", as_cmap=True), ax=ax[1])
ax[1].set_title("Average Growth Difference's\nInfluence of Services on Sectors")
del tps
# Layout
plt.tight_layout(pad=0)
plt.show()<jupyter_output><empty_output><jupyter_text>**How to Interpret: ** The color reflects the weight each service has on the sector. The weight is computed by averaging the growth of the sector by service, then normalizing these value by the sum of sectors growth. Each row sums to 1.**Interpretation: ** The general takeaway is that average growth suggests Risk and HR/Change are the major services, while growth rate is coming predominantely out of Technology while Risk plays a much smaller role.***Next I want to showcase my intial attempt at this viz and explain why it failed to simplify this large amount of information..<jupyter_code>f, ax = plt.subplots(1,2,figsize=[13,10], sharey=True)

# Stacked Bar Plot 1
tps = df.pivot_table(values=['mean'], index='Sector',columns='Service',aggfunc='mean')
tps.columns = tps.columns.droplevel()
tps = tps.div(tps.sum(1), axis=0)
tps.plot(kind='barh', stacked=True, ax=ax[0])
ax[0].legend(fontsize='large', loc='center left', bbox_to_anchor=(1, 0.5))
ax[0].set_xlabel("Average Growth Influence")
ax[0].set_title("Average Growth's\nInfluence of Services on Sectors")
ax[0].legend_.remove()
del tps

# Stacked Bar Plot 2
tps = diff.pivot_table(values=['mean'], index='Sector',columns='Service',aggfunc='mean')
tps.columns = tps.columns.droplevel()
tps = tps.div(tps.sum(1), axis=0)
tps.plot(kind='barh', stacked=True, ax=ax[1])
ax[1].legend(fontsize='large', loc='center left', bbox_to_anchor=(1, 0.5))
ax[1].set_xlabel("Average Growth Influence")
ax[1].set_title("Average Growth Diffrence's\nInfluence of Services on Sectors")
del tps
gc.collect()

plt.tight_layout(pad=0)
plt.show()<jupyter_output><empty_output><jupyter_text>**Criticism: ** While this viz attemps to communicate the same data as the heatmap above, the positioning on the bars make it difficult to compare influence between sectors and between metrics.**Deeper Dive into Service Influence on Sector over Time: ** Here is a resource to dive deeper into the over-time influence of services on each sector.<jupyter_code>plot_vars = sorted(melt_df["Sector"].unique())
f, ax = plt.subplots(len(plot_vars),2, figsize = [15,100], squeeze= True)

# By Specific Variable
for i,var in enumerate(plot_vars):
    # Column 1 Plots
    sns.pointplot(x= "Year", y= "Growth", hue= "Service", data= melt_df.loc[melt_df.Sector == var,:], ci= 25, ax= ax[i, 0])
    ax[i, 0].set_title("{}'s AvgGrowth by Year".format(var.title()))
    ax[i, 0].set_ylabel("Average Difference")
    ax[i, 0].legend_.remove()
    
    # Column 2 Plots
    sns.pointplot(x= "Year", y= "Growth", hue= "Service", data= diff_melt_df.loc[diff_melt_df.Sector == var,:], ci= 25, ax= ax[i, 1])
    ax[i,1].set_title("{}'s Avg Diff Growth by Year".format(var.title()))
    ax[i, 1].set_ylabel("Average Growth Difference")
    ax[i, 1].legend(fontsize='large', loc='center left', bbox_to_anchor=(1, 0.5))

# Layout
plt.tight_layout(pad=1)
plt.show()<jupyter_output><empty_output><jupyter_text>Q2: Which firm type is having the biggest influence on growth in the top 5 fastest growing countries between 2015-2019?<jupyter_code>top5_countries  = diff[["mean","Country"]].groupby("Country").mean().sort_values(by="mean", ascending= False)[:5]
print("Top 5 Fastest Growing Countries - Measured by Average Growth Difference by Year")
display(top5_countries.T)
print("Glancing at the African Continent figures..")
display(diff.loc[df.Region == "Africa", ["mean","Country"]].groupby("Country").mean().sort_values(by= "mean", ascending=False).T)<jupyter_output><empty_output><jupyter_text>Tremendous numbers coming out of Ethiopia. Since it grew 6 times faster than the runner up Canada, I am a little bit suspicious of the reliability of this source. Ideally I would follow up with data's source to ensure there is no mistake.<jupyter_code># Heatmap
f, ax = plt.subplots(1,1, figsize=[8,7], sharey=True)
tps = diff.loc[diff.Country.isin(top5_countries.index),:].pivot_table(values=['mean'], index='Country',columns='Firm Type',aggfunc='mean')
tps.columns = tps.columns.droplevel()
tps = tps.div(tps.sum(axis=1), axis=0)
tps = tps.loc[top5_countries.index]

# Plot Figure
sns.heatmap(tps, annot=True, fmt=".2f",cbar_kws={'label': 'Percent Influence'},cmap= sns.light_palette("red", as_cmap=True), ax=ax)
ax.set_title("Average Growth Difference's\nInfluence of Firm Type on Top 5 Countries")
ax.set_ylabel("Ranked Countries [Top to Bottom]")
del tps
# Layout
plt.yticks(rotation=0) 
plt.tight_layout(pad=0)
plt.show()<jupyter_output><empty_output><jupyter_text>**Take-Away: ** This heatmap suggests that across the board, Type P has the biggest influence on the top 5 growing countries. However, it is interesting to note that Type A dominates Canada, and Type M dominates Portugal.*** Q3: Can you forecast the current data set out to 2022? Use any methodology and external data sources you think are relevant.**Method:** Row-Wise Linear Models - Regression**Reasoning: **This is a very simple dataset with a limited number of variables. Furthermore, there is no seasonality, random-walk, or cyclical behavior to be found. Therefore, I believe that a simple linear will suffice to capture the upward trend. I hypothesize that using a polynomial term will improve the model by enabling it to capture the growth rate trend I have noticed. That being said, any significant disturbance in the system, such as a financial crisis, will render these projections obsolete, especially the polynomial projection, which is less robust than the linear prediction.<jupyter_code>df.head()
# Variable to Predict
y = df[['2013', '2014', '2015', '2016', '2017', '2018', '2019']]
years = len(y.columns)

# Scale years to low magnitude
X = np.array(range(1,8))

# Fit Line Function
from statistics import mean
def best_fit_slope_and_intercept(xs,ys):
    """
    Fit OLS Regression Line
    """
    m = (((mean(xs)*mean(ys)) - mean(xs*ys)) / ((mean(xs)*mean(xs)) - mean(xs*xs)))
    b = mean(ys) - m*mean(xs)
    return m, b
def predict_extrapolate(X, y):
    # Linear Regression from Scratch
    m, b = best_fit_slope_and_intercept(X,y)
    regression_line = [(m*x)+b for x in X]
    X_extrapolate = np.array(list(range(8, 11)))
    regression_line_future = [(m*x)+b for x in X_extrapolate]
    
    # Polynomial Regression with Numpy (2 Degrees)
    poly_coefs = np.polyfit(X,y, 2)
    ffit = np.polyval(poly_coefs, list(range(1,11)))

    # Linear Line
    plt.scatter(X,y,color='#003F72', label = "Training")
    plt.scatter(X_extrapolate,regression_line_future,color='g', label= "Linear Prediction")
    plt.plot(list(range(1,11)), regression_line + regression_line_future, color='g')
    
    # Polynomial Line
    # plt.scatter(X,y_train,color='#003F72', label = "Training")
    plt.plot(list(range(1,11)), ffit)
    plt.scatter(X_extrapolate,ffit[-3:],color='r', label= "2nd Polynomial Prediction")
    
    plt.legend(loc=4)
    plt.show()
for inter in range(10):
    random_index = random.randint(1,y.shape[0])
    print("Sample: ", random_index, "\n")
    dictionary = df.iloc[random_index,:7].to_dict()
    key = [x for x in dictionary]
    [print("{}: {}".format(k, dictionary[k])) for k in key]
    predict_extrapolate(X,np.array(y.iloc[random_index,:]))
    print("-------------------------------------------------\n")<jupyter_output><empty_output><jupyter_text>**Reflection on Polynomial Solution: ** From the sample case visualizations, the polynomial solution seems to capture the acceleration aspect of the trend, but I am concerned that it is extrapolating too aggresively. Perhaps a more ideal solution would take the average of the linear and polynomial solution in order to restrain the extrapolation, creating a more robust and conservative prediction.**Roll out Polynomial predictions to entire dataset..**<jupyter_code>def polynomial_fit_predict(X, y):
    poly_coefs = np.polyfit(X,y, 2)
    ffit = np.polyval(poly_coefs, list(range(8,11)))
    return ffit
%%time
df[["2020","2021","2022"]] = pd.DataFrame(df[['2013', '2014', '2015', '2016', '2017', '2018', '2019']].apply(
    lambda s: polynomial_fit_predict(X,s), axis=1).values.tolist())
print("Dataset with Prediction Sample:")
df.drop(["mean","stdev"],axis=1).sample(8)
print("Notebook Runtime: %0.2f Minutes"%((time.time() - notebookstart)/60))<jupyter_output><empty_output><jupyter_text>Manchester AirBnb Assignment_By <NAME>, January 2020_****Column Explanations****:**- **Property ID:** Airbnb Property ID- **Property Type:** Actual property type of the property being listed on Airbnb (Detached house, Flat etc)- **Listing Type:** The listing choice selected between Entire home, Apartment and Private room Bedrooms:** Number of bedrooms- **Reporting Month:** The month corresponding to the provided data- **Occupancy Rate:** Occupancy rate achieved in the given month- **Revenue (USD):** Revenue achieved in the given month (USD)- **Revenue (Native):** Revenue achieved in the given month (In currency of booking)- **ADR (USD):** Average daily revenue achieved in the given month (USD)- **ADR (Native):** Average daily revenue achieved in the given month (In currency of booking)- **Number of Reservations:** Number of reservations made for the listing in the given month- **Reservation Days:** Number of days listing was booked in the month- **Available Days:** Number of days listing was not booked in the month- **Blocked Days:** Number of days listing was blocked by the host in the month- **Active:** Whether property is active on Airbnb in the month or not- **Currency Native:** Currency of the booking being made- **Host ID:** Host identifier<jupyter_code>import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

import time
notebookstart = time.time()

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, cross_validate, cross_val_predict, KFold
from sklearn.metrics import mean_absolute_error
import lightgbm as lgb
from sklearn.preprocessing import StandardScaler

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import itertools

import math
import gc

%matplotlib inline
sns.set_style("whitegrid")

pd.options.display.max_rows = 999
pd.options.display.width = 300
pd.options.display.max_columns = 500

plot_title_size = 16
def plot_cloud(string, ax, title = "WordCloud", cmap = "plasma"):
    wordcloud = WordCloud(width=800, height=500,
                          collocations=True,
                          background_color="white",
                          colormap=cmap
                          ).generate(string)

    ax.imshow(wordcloud, interpolation='bilinear')
    ax.set_title(title,  fontsize=18)
    ax.axis('off')

def big_word_cloud(plot_df, plt_set, columns, figsize, cmap = "plasma"):
    """
    Iteratively Plot WordClouds
    """
    rows = math.ceil(len(plt_set)/columns)
    n_plots = rows*columns
    f,ax = plt.subplots(rows, columns, figsize = figsize)
    for i in range(0,n_plots):
        ax = plt.subplot(rows, columns, i+1)
        if i < len(plt_set):
            str_col = plt_set[i]
            string = " ".join(plot_df.loc[plot_df[str_col].notnull(),str_col].astype(str).str.lower().str.replace("none", "").str.title())
            string += 'EMPTY'
            ax = plt.subplot(rows, 2, i+1)
            plot_cloud(string, ax, title = "{} - {} Missing".format(str_col.title(), plot_df[str_col].isnull().sum()), cmap = cmap)
        else:
            ax.axis('off')
    plt.tight_layout(pad=0)

def big_count_plotter(plot_df, plt_set, columns, figsize, custom_palette = sns.color_palette("Dark2", 15), top_n = 15):
    """
    Iteratively Plot all categorical columns
    Has category pre-processing - remove whitespace, lower, title, and takes first 30 characters.
    """
    rows = math.ceil(len(plt_set)/columns)
    n_plots = rows*columns
    f,ax = plt.subplots(rows, columns, figsize = figsize)
    for i in range(0,n_plots):
        ax = plt.subplot(rows, columns, i+1)
        if i < len(plt_set):
            c_col = plt_set[i]
            plt_tmp = plot_df.loc[plot_df[c_col].notnull(),c_col].astype(str).str.lower().str.strip().str.title().apply(lambda x: x[:30])
            plot_order = plt_tmp.value_counts().index[:top_n]
            sns.countplot(y = plt_tmp, ax = ax, order = plot_order, palette = custom_palette)
            ax.set_title("{} - {} Missing".format(c_col.title(), plot_df[c_col].isnull().sum()), fontsize = plot_title_size)
            ax.set_ylabel("{} Categories".format(c_col.title()))
            ax.set_xlabel("Count")
        else:
            ax.axis('off')

    plt.tight_layout(pad=1)
    
def big_boxplotter(plot_df, plt_set, columns, figsize, custom_palette = sns.color_palette("Dark2", 15), quantile = .99):
    rows = math.ceil(len(plt_set)/columns)
    n_plots = rows*columns
    f,ax = plt.subplots(rows, columns, figsize = figsize)
    palette = itertools.cycle(custom_palette)
    for i in range(0,n_plots):
        ax = plt.subplot(rows, columns, i+1)
        if i < len(plt_set):
            cont_col = plt_set[i]
            plt_tmp = plot_df.loc[plot_df[cont_col].notnull(),cont_col].astype(float)
            col_max = plt_tmp.max()
            if quantile:
                plt_tmp = plt_tmp[plt_tmp < plt_tmp.quantile(quantile)]
            sns.boxplot(plt_tmp, color = next(palette))
            ax.set_title("{}\n{} Missing - {} Max".format(cont_col.title(), plot_df[cont_col].isnull().sum(), col_max),
                        fontsize = plot_title_size)
            ax.set_xlabel("Values (Upper Limit is {}th Quantile)".format(str(quantile)))
            ax.set_ylabel("{} Boxplot".format(cont_col))
        else:
            ax.axis('off')

    plt.tight_layout(pad=1)
    
def big_ts_count_plotter(plot_df, plt_set, columns, ts_var, figsize, custom_palette = sns.color_palette("Dark2", 15), top_n = 5):
    """
    Iteratively Plot all categorical columns over time
    """
    rows = math.ceil(len(plt_set)/columns)
    n_plots = rows*columns
    f,ax = plt.subplots(rows, columns, figsize = figsize)
    for i in range(0,n_plots):
        ax = plt.subplot(rows, columns, i+1)
        if i < len(plt_set):
            tmp_var = plt_set[i]
            plot_df.loc[plot_df[ts_var].notnull(),:].set_index(ts_var)['placeholder'].resample('1m').count().plot(ax=ax, label="Full Count")
            for tv in plot_df[tmp_var].dropna().value_counts().index[:top_n]:
                plot_df.loc[(plot_df[ts_var].notnull()) & (plot_df[tmp_var] == tv),:]\
                    .set_index(ts_var)['placeholder'].resample('1m').count().plot(ax=ax, label=tv)
            ax.legend(fontsize='large', loc='center left', bbox_to_anchor=(1, 0.5))
            ax.set_title("Monthly Count of {} Variable Through Time".format(tmp_var), fontsize = plot_title_size)
        else:
            ax.axis('off')

    plt.tight_layout(pad=1)

    
def big_ts_mean_plotter(plot_df, plt_set, columns, ts_var, std_multiplier, figsize, custom_palette = sns.color_palette("Dark2", 15)):
    """
    Iteratively Plot all continuous columns over time
    """
    rows = math.ceil(len(plt_set)/columns)
    n_plots = rows*columns
    f,ax = plt.subplots(rows, columns, figsize = figsize)
    palette = itertools.cycle(custom_palette)
    
    for i in range(0,n_plots):
        ax = plt.subplot(rows, columns, i+1)
        if i < len(plt_set):
            c = next(palette)
            tmp_var = plt_set[i]
            monthly_aggs = (plot_df
                           .set_index(ts_var)
                           .resample('1m')
                           .agg({tmp_var: ['mean','std']}))

            monthly_aggs.columns = pd.Index([e[0] +"_"+ e[1] for e in monthly_aggs.columns.tolist()])

            monthly_aggs = monthly_aggs.assign(
                upper = monthly_aggs[tmp_var + "_mean"] + (monthly_aggs[tmp_var + "_std"] * std_multiplier),
                lower = monthly_aggs[tmp_var + "_mean"] - (monthly_aggs[tmp_var + "_std"] * std_multiplier)
            )

            monthly_aggs[tmp_var+'_mean'].plot(ax=ax, color=c, label="{} Average".format(tmp_var))
            ax.fill_between(monthly_aggs.index, monthly_aggs.lower, monthly_aggs.upper, color=c, alpha=0.2)
            ax.legend(fontsize='large', loc='upper center', bbox_to_anchor=(0.5, -0.12))
            ax.set_title("Monthly Count of {} Variable Through Time".format(tmp_var), fontsize = plot_title_size)
        else:
            ax.axis('off')

    plt.tight_layout(pad=1)

def prepare_time_features(df, ts_var):
#     df['hour_of_day'] = df[ts_var].dt.hour
#     df['week'] = df[ts_var].dt.week
    df['month'] = df[ts_var].dt.month
    df["year"] = df[ts_var].dt.year
#     df['day_of_year'] = df[ts_var].dt.dayofyear
    df['week_of_year'] = df[ts_var].dt.weekofyear
#     df["weekday"] = df[ts_var].dt.weekday
    df["quarter"] = df[ts_var].dt.quarter
#     df["day_of_month"] = df[ts_var].dt.day
    
    return df

def time_slicer(plot_df, timeframes, value, agg_method, color="purple", figsize=[12,12]):
    """
    Function to count observation occurrence through different lenses of time.
    """
    f, ax = plt.subplots(len(timeframes), figsize = figsize)
    for i,x in enumerate(timeframes):
        plot_df.loc[:,[x,value]].groupby([x]).agg({value:agg_method}).plot(ax=ax[i],color=color)
        ax[i].set_ylabel(value.replace("_", " ").title())
        ax[i].set_title("{} by {}".format(value.replace("_", " ").title(), x.replace("_", " ").title()))
        ax[i].set_xlabel("")
        ax[i].set_ylim(0,)
    ax[len(timeframes)-1].set_xlabel("Time Frame")
    plt.tight_layout(pad=0)<jupyter_output><empty_output><jupyter_text>Load Data- Enforce data types- Create simple time features<jupyter_code>df = pd.read_csv("../input/housing-manchester/Airbnb_Manchester.csv")
print('DF shape: {} rows, {} columns'.format(*df.shape))

print("Number of unique Properties: {}".format(df['Property ID'].nunique()))
print("Number of unique Host: {}".format(df['Host ID'].nunique()))

display(df.sample(5))<jupyter_output><empty_output><jupyter_text>First impressions:- Panel Data, monthly time intervals. Between 2014 and 2019, which would very much reflect the dawn and boom of AirBnb.- Not much information about the actual housing products. - Operational features around revenue and occupancy.- GPS may be powerful feature to understand trendy neighbourhoods.TLDR - Fairly simple dataset with a decent number of observations (~200k) Data Cleaning and Exploration (40%)- Clean the dataset of unwanted data and outliers- Deal with missing data values Build Features:- How long Host / Property has been in action- Buckets for how enterpricing the Host is. (Owned properties)- Distance from Manchester City Center- Time features<jupyter_code>df['Reporting Month'] = pd.to_datetime(df['Reporting Month'])
df = prepare_time_features(df, ts_var = 'Reporting Month')
df['placeholder'] = 1
# Get the last active date of property
df = pd.merge(df,
         df.loc[df.Active == True,:].groupby('Property ID')['Reporting Month'].agg('max').rename("property_last_active").reset_index(),
         how='left', on = 'Property ID')
# Get the first active date of property
df = pd.merge(df,
         df.loc[df.Active == True,:].groupby('Property ID')['Reporting Month'].agg('min').rename("property_first_active").reset_index(),
         how='left', on = 'Property ID')
# Get the first active date of property
df = pd.merge(df,
         df.loc[df.Active == True,:].groupby('Host ID')['Reporting Month'].agg('max').rename("host_last_active").reset_index(),
         how='left', on = 'Host ID')
df = pd.merge(df,
         df.loc[df.Active == True,:].groupby('Host ID')['Reporting Month'].agg('min').rename("host_first_active").reset_index(),
         how='left', on = 'Host ID')

# Create features for Host and Properties around activity over time
df = df.assign(
    # Since First Active in days
    days_since_properties_first_active = (df['Reporting Month'] - df['property_first_active']).dt.days,
    days_since_host_first_active = (df['Reporting Month'] - df['host_first_active']).dt.days,
    # Until last Active in days
    days_until_properties_last_active = (df['property_last_active'] - df['Reporting Month']).dt.days,
    days_until_host_last_active = (df['host_last_active'] - df['Reporting Month']).dt.days,
    # Total life in days
    total_days_property_active = (df['property_last_active'] - df['property_first_active']).dt.days,
    total_days_host_active = (df['host_last_active'] - df['host_first_active']).dt.days,
    # Host Related Features
    host_owned_properties = df.groupby('Host ID')["Property ID"].transform('nunique')
)<jupyter_output><empty_output><jupyter_text>Filter Table- I have observed that around 40% of the dataset has non-active houses at a given month. These cases are void of revenue information, and do not contribute additional information. In fact, they provide a risk of misguiding the analysis.- The Country, State, City, Currency only have a single factor suiting Manchester. ZipCode is 100% Null. This will be ignored.- I will frame this analysis in terms of pounds, because this currency is local and will ignore exchange rates.<jupyter_code>print(df.Active.value_counts(normalize=True))
drop_cols =  [
    'Revenue (USD)',
    'ADR (USD)',
    'Country',
    'State',
    'City',
    'Zipcode',
    'Currency Native'
]

df = df.loc[df.Active == True,:].drop(drop_cols, axis=1)
print('DF shape: {} rows, {} columns'.format(*df.shape))
df.sample(5)<jupyter_output><empty_output><jupyter_text>Quick EDA<jupyter_code>sns.boxplot(df.groupby('Host ID')["Property ID"].agg('nunique'))
plt.title("Uniqued Properties owned by HostID")
plt.show()
count_plot_cols = [
    'Property Type',
    'Listing Type',
    'Bedrooms'
]

big_count_plotter(plot_df = df,
                  plt_set = count_plot_cols,
                  columns = 2,
                  figsize = [15,10])<jupyter_output><empty_output><jupyter_text>Main takeaways- House type is dominated by single bedroom Apartments or Houses.<jupyter_code>big_ts_count_plotter(plot_df = df.loc[df.Active == True,:],
                      plt_set = count_plot_cols,
                      columns = 2,
                      ts_var = 'Reporting Month',
                      figsize = [20,10])<jupyter_output><empty_output><jupyter_text>Main takeaways:- Appears to be a linear growth rate. This finding should be validated against unique properties.- Without the time-element, the top factors in categories appear more unequal than the reality. This is because the total difference is the sum of minor difference across (4 years * 12 months) ~48 montly periods<jupyter_code>f, ax = plt.subplots(1,2, figsize = [14,5])
tmp_plot = df.groupby('year').agg({"Host ID": "nunique", "Property ID": "nunique"})
tmp_plot.plot(ax=ax[0])
ax[0].set_ylim(0, )
ax[0].set_title("Number of Unique Host and Properties", fontsize = plot_title_size)

tmp_plot['Properties to Host Ratio'] = tmp_plot['Property ID'] / tmp_plot['Host ID']
tmp_plot['Properties to Host Ratio'].plot(ax=ax[1])
ax[1].set_ylim(1, 2)
ax[1].set_title("Properties to Host Ratio", fontsize = plot_title_size)
ax[1].set_ylabel("Properties to Host Ratio")
plt.tight_layout(pad=2)
plt.show()<jupyter_output><empty_output><jupyter_text>Main takeaways: - Over the last 4 years, there has been a minor increase in the average number of properties held by a host.- 2019 suggests a dip because the year is not complete yet.<jupyter_code>continuous_cols = [
    'Occupancy Rate',
    'Revenue (Native)',
    'ADR (Native)',
    'Number of Reservations',
    'Reservation Days',
    'Available Days',
    'Blocked Days',
    'days_since_properties_first_active',
    'days_since_host_first_active',
    'days_until_properties_last_active',
    'days_until_host_last_active',
    'total_days_property_active',
    'total_days_host_active'
]

big_boxplotter(plot_df = df,
               plt_set = continuous_cols,
               columns = 3,
               figsize = [20,25],
               quantile = .98)<jupyter_output><empty_output><jupyter_text>Main takeaways: - Host have been on the platform for an average of 2 years.- A property is available for 15 days on average a month, but also has a very wide variance.- A property is blocked for around 3 days on average a month.- Occupancy rate is around 35% on average- The revenue columns seem to have a outlier, perhaps a handful of mega-profitable houses. These should be double checked..<jupyter_code>df.loc[df['Revenue (Native)'] > 4000,:].sort_values(by='Revenue (Native)', ascending=False).iloc[:10]<jupyter_output><empty_output><jupyter_text>I just looks like there was an expensive Condominium that was rented out for 20 days at 7.8K pounds a day!<jupyter_code>continuous_cols = [
    'Occupancy Rate',
    'Revenue (Native)',
    'ADR (Native)',
    'Number of Reservations',
    'Reservation Days',
    'Available Days',
    'Blocked Days',
    'days_since_properties_first_active',
    'days_since_host_first_active',
    'days_until_properties_last_active',
    'days_until_host_last_active',
    'total_days_property_active',
    'total_days_host_active'
]

f, ax = plt.subplots(figsize = [8,5])
sns.heatmap(df.loc[:,continuous_cols].corr(method = 'spearman'),
            annot=False, fmt=".2f",cbar_kws={'label': 'Correlation Coefficient'},
            cmap="coolwarm",ax=ax, linewidths=.3, vmin=-1, vmax=1)
ax.set_title("Continuous Variables Correlation Matrix", fontsize = 14)
plt.show()<jupyter_output><empty_output><jupyter_text>Main takeaways: - Using Spearman rank for more robust (non-parametric) correlation test.- The revenue and availability features are highly correlated (Share data-generation process)- Days columns are highly correlated (Share data-generation process).- I'm surprised to see that the days property has been on the platform does not correlate with more occupancy. I would of thought that there would be a feedback loop where successful AirBnBs are less likely to drop off, and that successful AirBnbs get more comments and attention, hence more rents (network effect)<jupyter_code>upper_quantile = .90
lower_qunatile = .10

sns.jointplot(x="Longitude", y="Latitude", data=df.loc[
    (df['Longitude'] < df['Longitude'].quantile(upper_quantile)) &
    (df['Longitude'] > df['Longitude'].quantile(lower_qunatile)) &
    (df['Latitude'] < df['Latitude'].quantile(upper_quantile)) &
    (df['Latitude'] > df['Latitude'].quantile(lower_qunatile))
    ,["Longitude", "Latitude"]], kind="kde", color = 'green')
plt.title("2D Density Plot");<jupyter_output><empty_output><jupyter_text>Main takeaways: - There are two epicenter neighbourhoods- There are also some mild satelite areas in the general south west direction.<jupyter_code>continuous_cols = [
    'Occupancy Rate',
    'Revenue (Native)',
    'ADR (Native)',
    'Number of Reservations',
    'Reservation Days',
    'Available Days',
    'Blocked Days'
]

big_ts_mean_plotter(plot_df = df,
                      plt_set = continuous_cols,
                      columns = 2,
                      std_multiplier = .5,
                      ts_var = 'Reporting Month',
                      figsize = [20,25])<jupyter_output><empty_output><jupyter_text>Main takeaways: - There appears to be some seasonal patterns to the occupancy features. The pattern is an annual cycle of slow buildup until a peak in the fall, followed by a stabalizing drop.- Should followup to see whether there is an annual trend to observe.- There is a trend of aparments being available for a greater number of days over the years.<jupyter_code>time_slicer(plot_df = df,
            agg_method = 'mean',
            timeframes=['year', 'quarter', "month"],
            value = "Occupancy Rate", color="green", figsize = [10,7])<jupyter_output><empty_output><jupyter_text>Main takeaways: - Minor increase in occupancy rate over the years.- Annual cycle is lowest at .35 during new years, and .55 in september.<jupyter_code>time_slicer(plot_df = df,
            agg_method = 'sum',
            timeframes=['year', 'quarter', "month"],
            value = "Revenue (Native)", color="red", figsize = [10,7])<jupyter_output><empty_output><jupyter_text>Question of HotnessNow that I have explored the data, I am prepared to determine my dependent variable. My take on hotness is that it is primarily a reflection of what customers deem desirable. For this reason, many aspects of the revenue variables are disqualified, even though an economic argument could be made that the more money people are willing to spend, the more they must like a given property. This line of thinking may be distorted by the fact that ADR (Average Daily Revenue in given month) has an average of 50, but also a significant right tail that leads into the three hundred levels and over.This leaves the availability related variables. While occupancy rate could have high potential, it may also be distorted by apartments with high blocked days and average reservation days scenarios.This leaves me between reservation days and number of reservations. This is a tough trade off because on the one hand, reservation days suggests that available and rented apartments will climb to the top of the hotness metric, but it is at risk of scenarios where a single individuals rents an apartment over a long period of time. In this case, the motivations for renting are more likely tied to need for shelter.  Number of reservations, on the other hand, is representative of customers going through the apartment scoping process and selecting a given apartment X times in a given month. This product therefore most often persuades people to buy it. Data Aggregation (10%)- Aggregate the data by Property ID to create a table of yearly averages of the suitable fields (such as revenue, occupancy etc)- Convert the dataset that is in the form of monthly incomes to one which is yearly, and only keep data between 2018 and 2019<jupyter_code>columns_for_annual_sum = [
    'Revenue (Native)',
    'Number of Reservations',
    'Reservation Days',
    'Available Days',
    'Blocked Days'
]

static_columns = [
    'Property ID',
    'Property Type',
    'Listing Type',
    'Bedrooms',
    'Latitude',
    'Longitude',
    'year',
    'Active',
    'Host ID',
    'total_days_property_active',
    'total_days_host_active',
    'host_owned_properties']
model_df = df.loc[df.year == 2018,:].copy()
print('Model DF shape: {} rows, {} columns'.format(*model_df.shape))

agg_df = model_df.groupby(['year','Property ID']).agg({k:'sum' for k in columns_for_annual_sum}).reset_index()
print('agg_df DF shape: {} rows, {} columns'.format(*agg_df.shape))
year_df = pd.merge(
    agg_df,
    model_df.loc[:,static_columns].drop_duplicates(), on = ['year','Property ID'], how = 'left', validate = 'one_to_one')
print('year_df DF shape: {} rows, {} columns'.format(*year_df.shape))
display(year_df.sample(5))
del agg_df, model_df
# Rebuild lost features
year_df = year_df.assign(
    native_ADR = year_df['Revenue (Native)'] / year_df["Reservation Days"], # Revenue / Reservation days = ADR
    occupancy = year_df['Reservation Days'] / 365
).set_index('Property ID')

year_df['native_ADR'].fillna(0, inplace=True)
year_df['Reservation Days'].describe()<jupyter_output><empty_output><jupyter_text>Feature Engineering and Exploration (40%)- Normalise the data and carry out relevant transformations such as Integer encoding or One-Hot encoding if necessary- Engineer basic features (apart from the ones available) for predicting which properties are considered ‘hot’ (interpret this in your own way)<jupyter_code>dependent_variables  = 'Reservation Days'
categorical_variables = [
    'Property Type',
    'Listing Type'
]
continuous_variables = [
    'native_ADR',
    'host_owned_properties',
    'Bedrooms',
    'Latitude',
    'Longitude',
    'total_days_property_active',
    'total_days_host_active',
]
def dummy_variables(process_df, model_categorical):
    df = process_df.copy()
    dummy_df = pd.DataFrame()
    for cat_col in model_categorical:
        df[cat_col] = df[cat_col].str.lower().str.strip()
        print("{} Missing Values for {} - Dtype {}".format(df[cat_col].isnull().sum(), cat_col, df[cat_col].dtypes))
        df.loc[~df[cat_col].isin(df[cat_col].value_counts().index[:10]), cat_col] = cat_col + 'Other'
        
        
        dummies = pd.get_dummies(df[cat_col])
        dummies.columns = [x + "_" + cat_col for x in dummies.columns]
        dummy_df = pd.concat([dummy_df, dummies],axis = 1)
        
    return dummy_df

def standardscaler_forpd(process_df, model_continuous):
    df = process_df.loc[:,model_continuous].copy()
    for col in df.columns:
        scaler = StandardScaler()
        df[col] = scaler.fit_transform(df[col].values.reshape(-1, 1))
        
    return df
# Dummy variables for categorical..
dummy_df = dummy_variables(process_df = year_df, model_categorical = categorical_variables)
print(dummy_df.shape)

# Standard Scaler for Continuous..
cont_df = standardscaler_forpd(process_df = year_df, model_continuous = continuous_variables)
print(cont_df.shape)
display(cont_df.describe())

X = pd.concat([cont_df, dummy_df], axis =1)
y = year_df[dependent_variables].copy()
print("Dependent Variable Shape",y.shape)
print("Training Set Shape",X.shape)
del dummy_df, cont_df<jupyter_output><empty_output><jupyter_text>Modelling (10%)- Use these features plus the ones available in the data already and do some basic analysis and construct a Machine Learning model of your choice (whatever you think is suitable) for generating the prediction- Evaluate the model using appropriate evaluation metrics and explain your choice<jupyter_code>X['intercept'] = 1
model = LinearRegression(fit_intercept = False)

kf = KFold(n_splits=5, shuffle=True, random_state=42)
predicted = cross_validate(model, X, y, cv=kf, scoring = 'neg_mean_absolute_error')
print("CV MAE: {:.4f} +/- {:.4f}".format(abs(predicted['test_score'].mean()), predicted['test_score'].std()))

predicted = cross_val_predict(model, X, y, cv=kf)

# https://scikit-learn.org/stable/auto_examples/model_selection/plot_cv_predict.html#sphx-glr-auto-examples-model-selection-plot-cv-predict-py
fig, ax = plt.subplots()
ax.scatter(y, predicted, edgecolors=(0, 0, 0))
ax.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=4)
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title('Cross Validation Predictions\nAnalyze Model Residuels vs. Actual Out of Fold')
plt.show()<jupyter_output><empty_output><jupyter_text>Good to see low fluctuations across folds, however, looking at the prediction output, linear regression is unable to higher values of my dependent variable.<jupyter_code>print("Light Gradient Boosting Regressor: ")
lgbm_params =  {
    'task': 'train',
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': 'mae',
    "learning_rate": 0.05,
#     "num_leaves": 180,
#     "feature_fraction": 0.50,
#     "bagging_fraction": 0.50,
#     'bagging_freq': 4,
#     "max_depth": -1,
#     "reg_alpha": 0.3,
#     "reg_lambda": 0.1,
#     "min_split_gain":0.2,
#     "min_child_weight":10,
#     'zero_as_missing':True
                }
lgtrain = lgb.Dataset(X, y ,feature_name = "auto")

lgb_cv = lgb.cv(
    params = lgbm_params,
    train_set = lgtrain,
    num_boost_round=2000,
    stratified=False,
    nfold = 5,
    verbose_eval=50,
    seed = 23,
    early_stopping_rounds=75)
X_train, X_valid, y_train, y_valid = train_test_split(
    X, y, test_size=0.10, random_state=23)

lgtrain = lgb.Dataset(X_train, y_train)
lgvalid = lgb.Dataset(X_valid, y_valid)

lgb_reg = lgb.train(
    lgbm_params,
    lgtrain,
    num_boost_round=1000,
    valid_sets=[lgtrain, lgvalid],
    valid_names=['train','valid'],
    early_stopping_rounds=50,
    verbose_eval=100)
print("Model Evaluation Stage")
print('Validation Set MAE:', mean_absolute_error(y_valid, lgb_reg.predict(X_valid)))

# Feature Importance Plot
f, ax = plt.subplots(figsize=[7,10])
lgb.plot_importance(lgb_reg, max_num_features=50, ax=ax)
plt.title("Light GBM Feature Importance")
plt.show()<jupyter_output><empty_output><jupyter_text>LGBM is able to get a better score of 34 (CV score), which is 40% improvement over regression. It also provides a feature importance plot which helps to double check whether there are un-reasonable features in the mix (leaky).<jupyter_code># https://scikit-learn.org/stable/auto_examples/model_selection/plot_cv_predict.html#sphx-glr-auto-examples-model-selection-plot-cv-predict-py
fig, ax = plt.subplots()
ax.scatter(y_valid, lgb_reg.predict(X_valid), edgecolors=(0, 0, 0))
ax.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=4)
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title('Single Validation Predictions\nAnalyze Model Residuels')
plt.show()
sns.distplot(y_valid - lgb_reg.predict(X_valid))
plt.title("Distribution of Errors")
plt.xlabel("Residuals")
plt.show()
print("Notebook Runtime: %0.2f Minutes"%((time.time() - notebookstart)/60))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13283305720689992, 'fraction_numerical': 0.014878341941980964, 'mean_word_length': 4.718651309678618, 'pattern_counts': {'":': 20, '<': 193, '<?xml version=': 0, '>': 189, 'https://': 9, 'lorem ipsum': 0, 'www.': 7, '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': '45', 'dir': 'jupyter-structured-clean-dedup', 'id': '365691', 'license': 'MIT', 'n_tokens_mistral': 20963, 'n_tokens_neox': 19895, 'n_words': 5384, 'path': 'various-projects.ipynb', 'repo_name': 'nicapotato/KaggleNotebooks', 'size': '61581'}