<jupyter_start><jupyter_text>Sparkify - Capstone ProjectIn this notebook a capstone project for a Udacity Nanodegree is presented. Scope of WorkSparkify is a fictional music streaming application. Users can play songs, like/dislike songs, add them to playlist and add friends. There is a free and a paid user level. They can upgrade, downgrade or cancel their service at any time.The scope of this project is to do a churn prediction. Therefore users could get special offers before they quit the service.The given data set contains every interaction of the users in Sparkify. This notebook only uses a tiny subset (128MB) of the full data set available (12GB). It comprises data of 255 users starting at the 1th of October 2018 until the 3rd of December 2018. MethodPyspark is used to analyse the data and build three different ML models:* Logistic Regression* Naive Bayes* Random ForestFor the data visualization Matplotlib and Seaborn is used. Numpy and Pandas are imported for some data handling.<jupyter_code># install pyspark
!pip install pyspark
# connect to Google drive
from google.colab import drive
drive.mount("/content/gdrive")

gpath = "/content/gdrive/MyDrive/Colab Notebooks/Sparkify/"
# import libraries
from pyspark.sql import SparkSession
import pyspark.sql.functions as F
from pyspark.sql.types import IntegerType
from pyspark.sql import Window

from pyspark.ml import Pipeline
from pyspark.ml.feature import StringIndexer, VectorAssembler, StandardScaler
from pyspark.ml.classification import LogisticRegression, NaiveBayes, RandomForestClassifier
from pyspark.ml.tuning import CrossValidator, ParamGridBuilder
from pyspark.ml.evaluation import MulticlassClassificationEvaluator

import json

import datetime
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
# create a Spark session
spark = SparkSession.builder \
    .master("local[4]") \
    .appName("Spark Project") \
    .getOrCreate()
spark.sparkContext.getConf().getAll()<jupyter_output><empty_output><jupyter_text>1. Load and Clean DatasetIn this notebook, the mini-dataset file is `mini_sparkify_event_data.json`. The dataset is loaded and cleaned, checking for invalid or missing data in userids or sessionids.<jupyter_code># load json file and make it persist
df = spark.read.json(gpath + 'mini_sparkify_event_data.json')
df.persist()
# get size of df
print(df.count(), len(df.columns))<jupyter_output>286500 18<jupyter_text>Drop NaN for user Id<jupyter_code>df = df.dropna(how="any", subset=["userId"])
# get size of df
print(df.count(), len(df.columns))<jupyter_output>286500 18<jupyter_text>Drop duplicates<jupyter_code>df = df.dropDuplicates()
# get size of df
print(df.count(), len(df.columns))<jupyter_output>286500 18<jupyter_text>2. Exploratory Data Analysis (EDA)EDA is performed on the small subset of the data. Three new features for churn, downgrade phase and duration of each user  with Sparkify are created. Explore DataFirst of all basic properties like shape and data type of the dataset is explored. Then comparisons between churn and no-churn users are made to find significant differences.<jupyter_code># get first entry
df.toPandas().head(1).transpose()
# print data types of dataframe
df.printSchema()
# get statistics of numerical features
df.describe().toPandas()
# check for NaN
df.select([F.count(F.when(F.isnan(c), c)).alias(c) for c in df.columns]).toPandas()
# check for missing values in userId
df.where(df.userId == '').toPandas().head()
# check for missing values in sessionId
df.where(df.sessionId == '').toPandas().head()
# get unique 'page' entries
df.select('page').distinct().collect()<jupyter_output><empty_output><jupyter_text>2.1 Create 'churn' column where column=1 if page='Cancellation Confirmation'<jupyter_code># set a flag if 'page' = 'Cancellation Confirmation'
flag_cancellation_event = F.udf(lambda x: 1 if x == "Cancellation Confirmation" else 0, IntegerType())
df = df.withColumn("cancellation", flag_cancellation_event("page"))
# window over all timestamps in descending order
windowval = Window.partitionBy("userId").orderBy(F.desc("ts")).rangeBetween(Window.unboundedPreceding, 0)
# sum flags over window
df = df.withColumn("churn", F.sum("cancellation").over(windowval))<jupyter_output><empty_output><jupyter_text>2.2 Create Downgrade / Phase If user submits downgrade a flag is set. In the second step sum of flags is calculated using window function across descending timestamps.<jupyter_code># set a flag if 'page' = 'Submit Downgrade'
flag_downgrade_event = F.udf(lambda x: 1 if x == "Submit Downgrade" else 0, IntegerType())
df = df.withColumn("downgraded", flag_downgrade_event("page"))
# sum flags over window
df = df.withColumn("phase", F.sum("downgraded").over(windowval))<jupyter_output><empty_output><jupyter_text>2.3 Explore Data  Get number of churned / not churned users<jupyter_code># get distint user id's for churn and no-churn users
churn_abs = df.select('userId').where(df.churn==1).distinct().count()
not_churn_abs = df.select('userId').where(df.churn==0).distinct().count()
# create pie chart
labels = 'churn', 'no-churn'
sizes = [churn_abs, not_churn_abs]

fig1, ax1 = plt.subplots(figsize=(8, 6))
ax1.title.set_text('Fig1: Proportion of churn / no-churn users')
ax1.pie(sizes, labels=labels, autopct='%1.1f%%',
        shadow=True, startangle=90, colors=['tab:orange', 'tab:blue'])
ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.
plt.show()<jupyter_output><empty_output><jupyter_text>What's the average song playing length of churned / not churned users?<jupyter_code># aggregate average song playing length for both user groups
print("Churn users:")
df.where((df.churn==1) & (df.page=='NextSong')).select('length').describe().show()
print("No-churn users:")
df.where((df.churn==0) & (df.page=='NextSong')).select('length').describe().show()<jupyter_output>No-churn users:
+-------+-----------------+
|summary|           length|
+-------+-----------------+
|  count|           191714|
|   mean|249.2091353888089|
| stddev|98.62588300146908|
|    min|          0.78322|
|    max|       3024.66567|
+-------+-----------------+<jupyter_text>Did churned or not churned users downgrade from paid to free more often?<jupyter_code># aggregate average downgrade phase for both user groups
churn_avg_phase = df.where(df.churn==1).agg(F.avg(df.phase)).collect()[0][0]
no_churn_avg_phase = df.where(df.churn==0).agg(F.avg(df.phase)).collect()[0][0]
# create bar chart
x = [0, 1]
height = [no_churn_avg_phase, churn_avg_phase]
N = 2
ind = np.arange(N)    # the x locations for the groups
width = 0.35

fig2,ax1 = plt.subplots(figsize=(8, 6))
ax1.title.set_text('Fig2: Comparison of paid phases between churn / no-churn')
ax1.bar(x, height, width, color=['tab:blue', 'tab:orange'])
ax1.set_ylabel('Average paid phase [%]')
ax1.set_xticks(ind)
ax1.set_xticklabels(['no-churn', 'churn'])
plt.show()<jupyter_output><empty_output><jupyter_text>How often did churned / not churned users dislike a song or downgrade it on average (per user)?<jupyter_code># count 'Thumbs Down' + 'Downgrade' page for both user groups divided by each user quantity
print("Number of disliked songs of churn users:")
df.where((df.churn==1) & ((df.page=='Thumbs Down') | (df.page=='Downgrade'))).count() / churn_abs
print("Number of disliked songs of no-churn users:")
df.where((df.churn==0) & ((df.page=='Thumbs Down') | (df.page=='Downgrade'))).count() / not_churn_abs<jupyter_output>Number of disliked songs of no-churn users:<jupyter_text>How often did churned / not churned users liked a song or add a song to a playlist on average?<jupyter_code># count 'Thumbs Up' page for both user groups divided by each user quantity
print("Number of liked songs of churn users:")
df.where((df.churn==1) & ((df.page=='Thumbs Up') | (df.page=='Add To Playlist'))).count() / churn_abs
print("Number of liked songs of no-churn users:")
df.where((df.churn==0) & ((df.page=='Thumbs Up') | (df.page=='Add To Playlist'))).count() / not_churn_abs<jupyter_output>Number of liked songs of no-churn users:<jupyter_text>How many Songs did churned / not churned users play on average?<jupyter_code># count number of songs played by both user groups
song_count_df = df.where(df.page=='NextSong').groupby(['churn', 'userId']).count().toPandas()

# create boxplot
fig3,ax1 = plt.subplots(figsize=(8, 6))
sns.boxplot(x='churn', y='count', data=song_count_df)
ax1.set_title('Fig3: Avg Played Songs per churn / no-churn users')
ax1.set_xticklabels(['no-churn', 'churn'])
ax1.set_xlabel('Users')
plt.show()<jupyter_output><empty_output><jupyter_text>Distribution between 'paid' and 'free' users of churned users<jupyter_code># count 'Cancellation Confirmation' for both user groups
churn_free = df.where((df.page=='Cancellation Confirmation') & (df.level=='free')).count()
churn_paid = df.where((df.page=='Cancellation Confirmation') & (df.level=='paid')).count()

# create bar chart
labels = ['free', 'paid']
height = [churn_free, churn_paid]

fig4,ax1 = plt.subplots(figsize=(8, 6))

ax1.bar(x, height)
ax1.set_title("Fig4: Distribution of 'paid' / 'free' churned users")
ax1.set_ylabel('Churned Users')
ax1.set_xticks(ind)
ax1.set_xticklabels(['free', 'paid'])
plt.show()<jupyter_output><empty_output><jupyter_text>2.4 Calculating Statistics per Day<jupyter_code># get first timestamp in dataset:
min_ts = df.agg(F.min(df.ts)).collect()[0][0]
# get latest timestamp in dataset:
max_ts = df.agg(F.max(df.ts)).collect()[0][0]
# udf to get weekday of timestamp
get_weekday = F.udf(lambda x: datetime.datetime.fromtimestamp(x / 1000.0). weekday())
# create column for weekday
df = df.withColumn("wday", get_weekday(df.ts))
# head of data frame
df.head()<jupyter_output><empty_output><jupyter_text>How often played churned / not churned users songs during the week days on average?<jupyter_code># calculate average of churn user songs per weekday
churn_avg_df = df.where((df.churn==1) & (df.page == "NextSong")).groupby(df.wday)\
    .agg((F.count(df.wday) / churn_abs).alias('Songs per User')).orderBy(df.wday.cast("float")).collect()

churn_avg_values = []
for i in range(len(churn_avg_df)):
  churn_avg_values.append(churn_avg_df[i]['Songs per User'])

# calculate ration between weekdays and weekends
print(f"Ratio between week / weekend for churn users: {sum(churn_avg_values[:5]) / sum(churn_avg_values[-2:])}")
# calculate average of no-churn user songs per weekday
not_churn_avg_df = df.where((df.churn==0) & (df.page == "NextSong")).groupby(df.wday)\
    .agg((F.count(df.wday) / not_churn_abs).alias('Songs per User')).orderBy(df.wday.cast("float")).collect()

not_churn_avg_values = []
for i in range(len(not_churn_avg_df)):
  not_churn_avg_values.append(not_churn_avg_df[i]['Songs per User'])

# calculate ration between weekdays and weekends
print(f"Ratio between week / weekend for no-churn users: {sum(not_churn_avg_values[:5]) / sum(not_churn_avg_values[-2:])}")
# create stacked bar chart of avg songs per weekday for both groups
labels = ['Mo', 'Tue', 'Wed', 'Thu', 'Fr', 'Sa', 'Su']
width = 0.35

fig5, ax1 = plt.subplots(figsize=(8, 6))

ax1.bar(labels, not_churn_avg_values, width, label='no-churn')
ax1.bar(labels, churn_avg_values, width, bottom=not_churn_avg_values,
       label='churn')

ax1.set_ylabel('Avg Songs')
ax1.set_title('Fig5: Avg Played Songs per Weekday')
ax1.legend()

plt.show()<jupyter_output><empty_output><jupyter_text>Distribution of 'page' entries<jupyter_code># count all 'page' entries, sort descending
page_df = df.groupBy(df.page).count().sort(F.col("count").desc()).toPandas()
height = page_df['count'].tolist()
labels = page_df.page.tolist()
# create bar plot
fig6,ax1 = plt.subplots(figsize=(8, 6))
ax1.bar(labels, height)
ax1.set_title("Fig6: Distribution of 'page' entries")
plt.ylabel('Count [log scale]')
plt.xticks(rotation=90)
plt.yscale('log')
plt.show()
# head of data frame
page_df.head(10)<jupyter_output><empty_output><jupyter_text>2.5 Findings:* Groups of churn / not-churned users is unbalanced (23% / 77%).* Paid users are more presented then free users is both groups (Figure 2 and 4).* Not churned users play more songs on average than churn users.* Both groups play less songs on weekends. Ratio between week / weekend is greater for churn users (Figure 5).* Not-churned user have a longer usage time. Therefore all statistics about not churned users have a greater impact (see next chapter).  3. Feature EngineeringIn the first part of this section additional features are calculated to get a better presentation of the data. The main purpose of this work is to prepare user based data to train ML models. Therefore in the second to fourth parts new dataframes with one row per user id is generated using SQL queries. All features are transformed to numerical data.In the last part all subset are joined to the final dataset which will be used to train the models. Pearson coefficient is used to find correlated features and reduce them.  3.1 Create additional features  Calculate usage period per user in days'Usage period' is the difference between the first and the last timestamp for each user in the dataset.<jupyter_code># calculate usage period
df = df.withColumn("duration", F.round((F.col("ts").cast("long")/1000 - F.col("registration").cast("long")/1000)/86400, 4))<jupyter_output><empty_output><jupyter_text>average usage period in days for churned / not churned users<jupyter_code># aggregate average user duration for churn users
df.where(df.churn==1).agg(F.avg(df.duration)).show()
# aggregate average user duration for no-churn users
df.where(df.churn==0).agg(F.avg(df.duration)).show()<jupyter_output>+-----------------+
|    avg(duration)|
+-----------------+
|67.15457756826291|
+-----------------+<jupyter_text>truncate feature "userAgent"<jupyter_code># create temp view of data frame
df.createOrReplaceTempView("df")
# get one userAgent entry per user using sql
user_agent_short = spark.sql('''
                            SELECT DISTINCT userId, MAX(userAgent) userAgent
                            FROM df
                            GROUP BY userId
                            ''').toPandas()
# drop NaN from query
user_agent_short = user_agent_short.dropna()
def get_os(x):
  '''
  function to check which OS is used for each userId

  INPUT: 'userAgent' column

  OUTPUT: 1=Mac, 2=WIN, 3=Linux else 0

  '''
  if 'Mac' in x:
    return 1
  elif 'Windows' in x:
    return 2
  elif 'Linux' in x:
    return 3
  else:
    return 0

# create new column 'os' from 'userAgent' using get_os()
user_agent_short['os'] = user_agent_short['userAgent'].apply(get_os)
# drop 'userAgent'
user_agent_short = user_agent_short.drop('userAgent', axis=1)
# check for 0
user_agent_short[user_agent_short.os == 0]
# create pyspark data frame
user_agent_py = spark.createDataFrame(user_agent_short)
# show head of pyspark data frame
user_agent_py.groupby(user_agent_py.os).count().show()<jupyter_output>+---+-----+
| os|count|
+---+-----+
|  1|  102|
|  3|   12|
|  2|  111|
+---+-----+<jupyter_text>3.2 Create subset dataframe of relevant features<jupyter_code># select relevant features from df
df_feat = df.select(['artist', 'gender', 'itemInSession', 'length', 'level', 'location', 'method', 'page',\
                     'registration', 'sessionId', 'song', 'status', 'ts', 'userId', 'churn', 'phase', 'wday',\
                     'duration'])
# convert userId, wday to integer
df_feat = df_feat.withColumn("userId", df_feat.userId.cast(IntegerType()))
df_feat = df_feat.withColumn("wday", df_feat.wday.cast(IntegerType()))
# head of data frame
df_feat.head()<jupyter_output><empty_output><jupyter_text>3.3 Convert string features to numerical features  Convert 'gender' column to integer<jupyter_code># udf to convert gender feature to numerical data type
convert_gender = F.udf(lambda x: 1 if x == 'M' else 0, IntegerType())
# create new column 'gender_num' using udf
df_feat = df_feat.withColumn('gender_num', convert_gender(df_feat.gender))<jupyter_output><empty_output><jupyter_text>'length': Replace None with 0<jupyter_code>df_feat = df_feat.fillna(0, subset=['length'])<jupyter_output><empty_output><jupyter_text>Convert 'level' column to integer, get lastest level<jupyter_code># udf to convert level feature to numerical data type
convert_level = F.udf(lambda x: 1 if x == 'paid' else 0, IntegerType())
# create new column 'level_num' using udf
df_feat = df_feat.withColumn('level_num', convert_level(df_feat.level))
# get lastest level of each user
levels = df_feat.select(['userId', 'level_num', 'ts']).orderBy(F.desc('ts')).dropDuplicates(['userId'])\
    .select(['userId', 'level_num'])<jupyter_output><empty_output><jupyter_text>Convert 'location' column to Divison Columnssource: https://github.com/cphalpert/census-regions/blob/master/us%20census%20bureau%20regions%20and%20divisions.csv<jupyter_code>df_feat.select(['location']).head(5)
# import csv file
region = spark.read.csv(gpath + 'devisions.csv', header=True)
# udf to get state code of location string if string is not None
find_loc = F.udf(lambda x: x.split(', ')[1].split('-')[0] if x is not None else x)
# udf to lower state code and create new column name
loc_name = F.udf(lambda x: 'state_'+x.replace(' ', '_').lower())

# create new data frame with unique state codes of df_feat
locations = df_feat.withColumn('location', find_loc(df_feat.location))\
        .select(['userId', 'location']).filter(df_feat.location.isNotNull()).dropDuplicates(subset=['userId'])

# join division name from region csv file
locations = locations.join(region, locations.location == region['State Code'], how='left')\
    .select(['userId', F.col("Division").alias("location")])\
    .withColumn('location', loc_name('location'))

locations = locations.groupBy('userId').pivot('location').agg(F.count('location')).fillna(0)
# show locations data frame
print(f"{locations.count()} rows.")
locations.toPandas().head(5)<jupyter_output>225 rows.<jupyter_text>Convert 'method' column to integer<jupyter_code>df_feat.select('method').distinct().show()
# udf to convert 'method' column to numerical data type
convert_method = F.udf(lambda x: 1 if x == 'PUT' else 0, IntegerType())
# create new column 'method_num' using udf
df_feat = df_feat.withColumn('method_num', convert_method(df_feat.method))<jupyter_output><empty_output><jupyter_text>3.4 Create Datasets with one row per user<jupyter_code># create temp view of df_feat
df_feat.createOrReplaceTempView("df_feat")<jupyter_output><empty_output><jupyter_text>User Dataframe<jupyter_code># create a user dataframe using sql query
user_num = spark.sql(
    '''
    SELECT userId, AVG(gender_num) gender, AVG(length) length,
    AVG(itemInSession) itemInSession , AVG(method_num) method, MIN(registration) registration, AVG(status) status,
    MAX(churn) churn, AVG(wday) avg_wday, MAX(duration) user_duration
    FROM df_feat
    GROUP BY userId
    '''
)

print(f"{user_num.count()} rows.")
user_num.toPandas().head(5)
# check for NaN
user_num.select([F.count(F.when(F.isnan(c), c)).alias(c) for c in user_num.columns]).toPandas()<jupyter_output><empty_output><jupyter_text>Pages Dataframe<jupyter_code># create a dataframe containing counts of all pages per userId
page_name = F.udf(lambda x: 'page_' + x.replace(" ", "_").lower())

user_pages = df_feat.select(['userId', 'page']).withColumn('page', page_name(df_feat.page))\
    .groupBy(['userId']).pivot("page").agg(F.count('page'))

# drop userId == 0 (NaN), fill rest of NaN with 0
user_pages = user_pages.where(F.col('userId').isNotNull()).fillna(0)

print(f"{user_pages.count()} rows.")
user_pages.toPandas().head(5)
# calculate average to check if one column only contains Null
user_pages.select(*[F.mean(c).alias(c) for c in user_pages.columns]).show()
# drop columns with average of 0
cols = ('page_login', 'page_register', 'page_submit_registration')
user_pages = user_pages.drop(*cols)
# drop 'Cancellation Confirmation' and 'Cancel'
# they are equal to churn, which would be a data leakage for the models
user_pages = user_pages.drop('page_cancellation_confirmation').drop('page_cancel')
# check for NaN
user_pages.select([F.count(F.when(F.isnan(c), c)).alias(c) for c in user_pages.columns]).toPandas()<jupyter_output><empty_output><jupyter_text>Artists / Song Dataframe<jupyter_code># count number of distinct artists, number of songs played and number of distinct songs
artists_num = spark.sql(
    '''
    SELECT userId, COUNT(DISTINCT artist) num_artists_dist, 
    COUNT(song) num_songs, COUNT(DISTINCT song) num_songs_dist
    FROM df_feat
    GROUP BY userId
    '''
)

print(f"{artists_num.count()} rows.")
artists_num.toPandas().head(5)
# check for NaN
artists_num.select([F.count(F.when(F.isnan(c), c)).alias(c) for c in artists_num.columns]).toPandas()<jupyter_output><empty_output><jupyter_text>3.5 Joining all features / reduce features<jupyter_code># join all data frame to get the final dataset of user data
users_dataset = user_num.join(locations, ['userId']).join(levels, ['userId'])\
            .join(user_pages, ['userId']).join(artists_num, ['userId']).join(user_agent_py, ['userId'])

users_dataset.toPandas().head(5)
# get size
print(users_dataset.count(), len(users_dataset.columns))<jupyter_output>225 41<jupyter_text>Calculate Pearson coefficient, drop features with p > 0.8To reduce the number of features Pearson coefficient is used to calculuate correlation between all features and the churn property.<jupyter_code># convert dataset to Pandas data frame
feat_pd = users_dataset.toPandas()
# calculate all Pearson coefficients
corr_matrix = feat_pd.corr()
# head of correlation matrix
corr_matrix.head()
# get lower triangle of corr matrix without diagonal
lower_tri = corr_matrix.where(np.tril(np.ones(corr_matrix.shape),k=-1).astype(np.bool))
# create list of columns (=features) which are > 0.8
to_drop = [col for col in lower_tri.columns if any(np.abs(lower_tri[col]) > 0.8)]
print(to_drop)<jupyter_output>['length', 'registration', 'page_add_friend', 'page_add_to_playlist', 'page_downgrade', 'page_help', 'page_home', 'page_logout', 'page_nextsong', 'page_settings', 'page_thumbs_down', 'page_thumbs_up', 'num_artists_dist', 'num_songs']<jupyter_text>Get Columns to Keep<jupyter_code># create list which features which are not dropped
cols_to_keep = feat_pd.columns.drop(to_drop).tolist()
# reduce correlation matrix containing only cols_to_keep
corr_matrix_remain = feat_pd[cols_to_keep].corr()
# get shape of corr matrix
corr_matrix_remain.shape<jupyter_output><empty_output><jupyter_text>Plot correlation coefficients<jupyter_code># Generate a mask for the lower triangle
mask = np.zeros_like(corr_matrix_remain, dtype=bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 11))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(230, 20, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr_matrix_remain, mask=mask, cmap=cmap, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})
ax.set_title("Fig7: Diagonal correlation matrix")
plt.show()
# keep only features with low correlation + churn
# reduce dataset to cols_to_keep
dataset = users_dataset.select(cols_to_keep)
# check for missing values / get shape / show columns
dataset_pd = dataset.toPandas()
missing_count = dataset_pd.isnull().sum().sum()
print(f"Missing values in dataset: {missing_count}.")
print(f"Shape of dataset: {dataset_pd.shape}")
print(f"Columns of dataset: {dataset_pd.columns}")<jupyter_output>Missing values in dataset: 0.
Shape of dataset: (225, 27)
Columns of dataset: Index(['userId', 'gender', 'itemInSession', 'method', 'status', 'churn',
       'avg_wday', 'user_duration', 'state_east_north_central',
       'state_east_south_central', 'state_middle_atlantic', 'state_mountain',
       'state_new_england', 'state_pacific', 'state_south_atlantic',
       'state_west_north_central', 'state_west_south_central', 'level_num',
       'page_about', 'page_error', 'page_roll_advert', 'page_save_settings',
       'page_submit_downgrade', 'page_submit_upgrade', 'page_upgrade',
       'num_songs_dist', 'os'],
      dtype='object')<jupyter_text>4. ModelingThe data set is split into train, test and validation sets. Each model is trained and tested. The best model will be validated using the validation set. The following models are used:**Logistic Regression**It is used for binary variables. It learns a linear relationship from the given data and then uses the Sigmoid function to implement a non-linearity. It is easy to implement, interpret and very efficient to train.[https://www.geeksforgeeks.org/advantages-and-disadvantages-of-logistic-regression/]**Naive Bayes**This model uses the Bayes' Theorem and assumes that all features are independent from each other thus there is no correlation between them. This assumption makes it very fast. It also works well with data with many features.[https://towardsdatascience.com/naive-bayes-classifier-explained-50f9723571ed]**Random Forest Classifier**This model builds many randomly generated dicision trees and averages over all predictions. It is relatively slow in training but has a lower risk of overfitting. This algorithm tends to be biased when dealing with categorical features.[https://towardsai.net/p/machine-learning/why-choose-random-forest-and-not-decision-trees]<jupyter_code># rename churn column to 'label', convert to float
dataset = dataset.withColumn('label', dataset['churn'].cast('float')).drop('churn')

# get feature cols
feature_cols = dataset.drop('label').drop('userId').columns
# split into train, test, validation sets
train, test, validation = dataset.drop('userId').randomSplit([0.6, 0.2, 0.2], seed=42)<jupyter_output><empty_output><jupyter_text>Evaluation F1-score is choosen for the evaluation metric. Because the data is unbalanced (churn: 23%, not churn: 77%). Therefore accuracy would be relatively high, if the model would always predict "not churn". F1-score is a trade off between precision (positive prediction) and recall (sensitivity).$$ F_{1} = \frac{2 * precision * recall}{precision + recall} $$  4.1 Logistic Regression Feature ScalingIt is recommended for regression to scale the features so that the predictors have a mean of zero. This makes the intercept term easier to interpret as the expected value of Y when the predictor values are set to their means.https://www.atoti.io/articles/when-to-perform-a-feature-scaling/:~:text=Feature%20scaling%20is%20a%20method,during%20the%20data%20preprocessing%20step.<jupyter_code># build pipeline
assembler = VectorAssembler(inputCols=feature_cols, outputCol="features")
scaler = StandardScaler(inputCol="features", outputCol="scaledFeatures", withStd=True, withMean=True)
lr = LogisticRegression(labelCol="label", featuresCol="scaledFeatures", maxIter=10, regParam=0.0, elasticNetParam=0)
pipeline = Pipeline(stages=[assembler, scaler, lr])
# param Grid
paramGrid = ParamGridBuilder().build()

# Cross validation
crossval = CrossValidator(
    estimator=pipeline,
    estimatorParamMaps=paramGrid,
    evaluator=MulticlassClassificationEvaluator(metricName='f1'),
    numFolds=3
)
# train model
cvLrModel = crossval.fit(train)
# evaluate on test data
pred_Lr = cvLrModel.transform(test)
evaluator = MulticlassClassificationEvaluator(predictionCol="prediction", labelCol="label")

accu = pred_Lr.where(pred_Lr.prediction == pred_Lr.label).count() / pred_Lr.count()
print(f"Accuracy: {round(accu, 4)}")
print("F1-score: ", evaluator.evaluate(pred_Lr, {evaluator.metricName: "f1"}))<jupyter_output>Accuracy: 0.7963
F1-score:  0.7906013540422143<jupyter_text>4.2 Naive Bayes Feature ScalingThis algorithm is not distance based, but probability based. Therefore it does not require feature scaling.https://medium.datadriveninvestor.com/naive-bayes-classifier-essential-things-to-know-c0d7d30b2954<jupyter_code># build pipeline
assembler = VectorAssembler(inputCols=feature_cols, outputCol="features")
Nb = NaiveBayes(labelCol="label", featuresCol="features", smoothing=1.0)
pipeline = Pipeline(stages=[assembler, Nb])
# param Grid
paramGrid = ParamGridBuilder().build()

# Cross validation
crossval = CrossValidator(
    estimator=pipeline,
    estimatorParamMaps=paramGrid,
    evaluator=MulticlassClassificationEvaluator(metricName='f1'),
    numFolds=3
)
# train model
cvNbModel = crossval.fit(train)
# evaluate on test data
pred_Nb = cvNbModel.transform(test)#.select('label', 'prediction')
evaluator = MulticlassClassificationEvaluator(predictionCol="prediction", labelCol="label", metricName='f1')

accu = pred_Nb.where(pred_Nb.prediction == pred_Nb.label).count() / pred_Nb.count()
print(f"Accuracy: {round(accu, 4)}")
print("F1-score: ", evaluator.evaluate(pred_Nb, {evaluator.metricName: "f1"}))<jupyter_output>Accuracy: 0.4815
F1-score:  0.5455337690631807<jupyter_text>4.3 Random ForestThe Random Forest Classifier is a decision tree based algorithm which is not affected by feature scaling. Therefore no scaler is used here.<jupyter_code># build pipeline
assembler = VectorAssembler(inputCols=feature_cols, outputCol="features")
Rf = RandomForestClassifier(labelCol="label", featuresCol="features", maxDepth=5, numTrees=30, seed=42)
pipeline = Pipeline(stages=[assembler, Rf])
# param Grid
paramGrid = ParamGridBuilder().build()    

# Cross validation
crossval = CrossValidator(
    estimator=pipeline,
    estimatorParamMaps=paramGrid,
    evaluator=MulticlassClassificationEvaluator(metricName='f1'),
    numFolds=3
)
# train model
cvRfModel = crossval.fit(train)
# evaluate on test data
pred_Rf = cvRfModel.transform(test)
evaluator = MulticlassClassificationEvaluator(predictionCol="prediction", labelCol="label")

accu = pred_Rf.where(pred_Rf.prediction == pred_Rf.label).count() / pred_Rf.count()
print(f"Accuracy: {round(accu, 4)}")
print("F1-score: ", evaluator.evaluate(pred_Rf, {evaluator.metricName: "f1"}))<jupyter_output>Accuracy: 0.8889
F1-score:  0.8589569160997731<jupyter_text>Hyperparamter TuningLogistic Regression implies a linear relationship of the features to the classes. Naive Bayes implies no correlation between all features. Both assumptions do not apply to this case. Therefore the scores show not the best results.The Random Forest Classifier has the highest F1-score. Therefore the hyperparameter of this model are tuned using the parameter grid.* Max. Depth: Overfitting should be avoided so I only increase the depth in steps of 2.* Number of Trees: The trees are randomly generated. Thus overfitting is not a risk here.The best model is then tested on the validation data.<jupyter_code># param Grid
paramGrid = ParamGridBuilder()\
    .addGrid(Rf.maxDepth,[5, 7]) \
    .addGrid(Rf.numTrees,[30, 50, 70]) \
    .build()
    

# Cross validation
crossval = CrossValidator(
    estimator=pipeline,
    estimatorParamMaps=paramGrid,
    evaluator=MulticlassClassificationEvaluator(metricName='f1'),
    numFolds=3
)
# train model
cvRfModel = crossval.fit(train)
# average metrics of each parameter setting on the training data:
cvRfModel.avgMetrics
# evaluate on validation data
pred_Rf = cvRfModel.transform(validation)
evaluator = MulticlassClassificationEvaluator(predictionCol="prediction", labelCol="label")

accu = pred_Rf.where(pred_Rf.prediction == pred_Rf.label).count() / pred_Rf.count()
print(f"Accuracy: {round(accu, 4)}")
print("F1-score: ", evaluator.evaluate(pred_Rf, {evaluator.metricName: "f1"}))
# get parameters of best model
bestModel = cvRfModel.bestModel
print("Max. depth of best model: ", bestModel.stages[-1]._java_obj.getMaxDepth())
print("Num. of trees of best model: ", bestModel.stages[-1]._java_obj.getNumTrees())<jupyter_output>Max. depth of best model:  5
Num. of trees of best model:  50<jupyter_text>4.4 Feature Importance<jupyter_code># get estimation of feature importances of best model
feature_imp = bestModel.stages[-1].featureImportances

# create pd dataframe of feature importances with labels
featImp_df = pd.DataFrame(data={'label': feature_cols, 'imp': feature_imp.toArray().tolist()})

# sort in ascending order
featImp_df.sort_values('imp', inplace=True)
# create bar chart
labels = featImp_df.label
width = featImp_df.imp.values
y = range(len(feature_cols))
ind = np.arange(len(feature_cols))
height = 0.8

my_cmap = plt.get_cmap('Wistia')
rescale = lambda x: (x - np.min(x)) / (np.max(x) - np.min(x))

fig8,ax1 = plt.subplots(figsize=(12, 8))

ax1.barh(y, width, height, color=my_cmap(rescale(width)))
ax1.set_title("Fig8: Feature Importances of Random Forest Model")
ax1.set_yticks(ind)
ax1.set_yticklabels(labels)
plt.show()
df_feat.select(['page', 'method']).distinct().toPandas().sort_values('method').head(30)
# investigate average method of churn and no-churn users
print("Average method of churn/no-churn users: ")
users_dataset.groupby(['churn']).agg(F.avg('method')).toPandas().head()
# investigate average weekday of churn and no-churn users
print("Average weekday of churn/no-churn users: ")
users_dataset.groupby(['churn']).agg(F.avg('avg_wday')).toPandas().head()
# investigate number of roll adverts in both groups
print("Average number of roll advertisment of both groups: ")
users_dataset.groupby(['churn']).agg(F.avg('page_roll_advert')).toPandas().head()
# investigate number of churn users in Middle Atlantic
print('Churn users in Middle Atlantic in %:')
round(users_dataset.where((users_dataset.state_middle_atlantic == 1) & (users_dataset.churn == 1)).count() \
  / users_dataset.where(users_dataset.state_middle_atlantic == 1).count() * 100, 2)<jupyter_output>Churn users in Middle Atlantic in %:<jupyter_text>5. ConclusionUser data from Sparkify was used to train a Random Forest Classifier to predict churn users. F1-score metric was used to get high precision as well as recall.The final model has a acceptable F1-score of about 81.75% on the validation data. The model was only trained on a small subset. Using a spark cluster and the whole 12GB data set would probably increase the score significantly.This model could be trained to predict the churn users on a weekly base and then send these users discounts or personalized offers. The impact of these actions could then be analysed with an A/B test.<jupyter_code># convert notebook to html file
#%%shell
#jupyter nbconvert --to html /content/Sparkify.ipynb<jupyter_output>[NbConvertApp] Converting notebook /content/Sparkify.ipynb to html
[NbConvertApp] Writing 781451 bytes to /content/Sparkify.html
{'fraction_non_alphanumeric': 0.1274424537600742, 'fraction_numerical': 0.01635066968168377, 'mean_word_length': 5.749168460183917, 'pattern_counts': {'":': 1, '<': 147, '<?xml version=': 0, '>': 150, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, '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': '41', 'dir': 'jupyter-structured-clean-dedup', 'id': '633215', 'license': 'W3C', 'n_tokens_mistral': 11690, 'n_tokens_neox': 11200, 'n_words': 3138, 'path': 'Sparkify.ipynb', 'repo_name': 'LollaPie/Sparkify', 'size': '34494'}