Data science project for ANLY-501 | A data science tale about a coffee company.

This notebook provides a project report for the data science project done as part of ANLY-501. The data science project intends to show all stages of the data science pipeline. This notebook facilitates that by organizing different phases into different sections and having the code and visualizations all included in one place. Install Jupyter notebook software available from http://jupyter.org/ to run this notebook. All code for this project is available on Github as https://github.com/aarora79/sb_study and there is also a simple website associated with it https://aarora79.github.io/sb_study/. The entire code (including the generated output) can be downloaded as a zip file from Github via this URL https://github.com/aarora79/sb_study/archive/master.zip. To run the code, simply unzip the master.zip and say "python main.py". The code requires that Python 3.5.2|Anaconda 4.1.1 be installed on the machine. Also, since the code tries to get the datasets in realtime so an Internet connection is also required.

[Updated for Project3 11/30/2016, starting from section Predictions based on classification].

Data Science Problem

Starbucks Corporation is an American coffee company and coffeehouse chain with more than 24,000 stores across the world. This projects intends to explore the following data science problems:

1. Exploratory Data Analysis (EDA) about Starbucks store locations for example geographical distribution of stores by country, region, ownership model, brand name etc.
2. Find a relationship between Starbucks data with various economic and human development indices such as GDP, ease of doing business, rural to urban population ratio, literacy rate, revenue from tourist inflow and so on.
3. Predict which countries where Starbucks does not have a store today are most suitable for having Starbucks stores (in other words in which country where Starbucks does not have a presence should Starbucks open its next store and how many).

This problem is important because it attempts to provide a model (using Starbucks as an example) which can be applied to any similar business (say in the food and hospitality industry) to predict where it could expand globally. It provides insight about where the business is located currently by bringing out information like say the total number of stores (48) found in the entire continent of Africa is way less than number of stores (184) found just on U.S. Airports.

Potential Analysis that Can Be Conducted Using Collected Data

The data used as part of this project is obtained from two sources.

1. Starbucks store location data is available from https://opendata.socrata.com via an API. The Socrate Open Data API (SODA) end point for this data is https://opendata.socrata.com/resource/xy4y-c4mk.json

2. The economic and urban development data for various countries is available from the World Bank(WB) website. WB APIs are described here https://datahelpdesk.worldbank.org/knowledgebase/topics/125589-developer-information. The API end point for all the World Development Indicators (WDI) is http://api.worldbank.org/indicators?format=json. Some examples of indicators being collected include GDP, urban to rural population ratio, % of employed youth (15-24 years), international tourism receipts (as % of total exports), ease of doing business and so on.

The possible directions/hypothesis based on collected data (this is not the complete list, would be expanded as the project goes on):

1. EDA about Starbucks store locations (for example):

   • What percentage of stores exists in high income, high literacy, high urban to rural population ratio European countries Vs say high population, rising GDP, low urban to rural population Asian countries.
   • Distribution of stores across geographies based on type of ownership (franchisee, joint venture etc.), brand name etc.
   • Which country, which city has the most Starbucks stores per 1000 people.
   • Is there a Starbucks always open at any UTC time during a 24hour period i.e. you can always find some Starbucks store open time at any given time somewhere in some timezone around the world.

2. Data visualization of the Starbucks store data:

   • World map showing starbucks locations around the world.
   • Heat map of the world based on the number of Starbucks store in a country.
   • Frequency distribution of Starbucks store by city in a given country. Does this distribution resemble any wellll known statistical distribution.
   • Parallel coordinates based visualization for number of stores combined with economic and urban development indicators.

3. Machine learning model for predicting number of Starbucks store based on various economic and urban development indicators. This could then be used to predict which countries where Starbucks does not have a presence today would be best suited as new markets. For example, model the number of Starbucks location in a country based on a) ease of doing business, b) GDP, c) urban to rural population, d) employment rate between people in the 15-24year age group, e) type of government, f) access to 24hour electricity and so on and so forth.

Data Issues

The data used for this project is being obtained via APIs from the Socrata web site and the World Bank website and is therefore expected to be relatively error free (for example as compared to the same data being obtained by scraping these websites). Even so, the data is checked for quality and appropriate error handling or even alternate mechanisms are put in place to handle errors/issues with the data.

The subsequent sections of this notebook provide a view of the data for the two datasets used in this project and also provide the data quality scores (DQS) for both the datasets.

Datasets (Starbucks and WorldBank WDI)

The Python code for the SB Study (SBS) is run offline and the results (along with the code) are periodically pushed to Github. The code here simply downloads the files from Github to show the results and how the data looks like.

Data Quality Scores (DQS)

To check for dataquality two types of checks were done on both the datasets.

1. Missing data: any cell in the dataset which was empty was counted as missing data. This is checked for all features in both the datasets.
2. Invalid data: validation checks were done on many (9 fields in the Starbucks dataset, and all 83 fields in the Worldbank dataset). These checks were both generic (validation for numeric data, absence of special characters, validation of timestamp data, latitude and longitude validation, timezone offset and so on and so forth) as well context specific (store-id has to be unique).

Two metrices are derived to quantify missing data and invalid data. These are raw score and adjusted score.

1. Raw Score for Missing Data: this is the percentage of data that is not missing, simply speaking this is the count of non-empty cells Vs total number of cells expressed as a percentage. The higher this score the cleaner the dataset is from the perspective of missing data. This score is available for both Starbucks and WorldBank datasets.

2. Adjusted Raw Score for Missing Data: this is the percentage of data that is not missing for Mandatory Features (i.e. features without which the record would have to be discarded), simply speaking this is the count of non-empty cells in mandatory columns Vs total number of cells expressed as a percentage. In case all features are mandatory then the adjusted raw score and the raw score are the same. The higher this score the cleaner the dataset is from the perspective of missing data. This score is available for both Starbucks and WorldBank datasets.

3. Raw Score for Invalid Data: this is the percentage of data that is not invalid, simply speaking this is the count of cells containing valid data Vs total number of cells expressed as a percentage. The higher this score the cleaner the dataset is from the perspective of invalid data. Validity checks are both generic and context specific as defined above. This score is available for both Starbucks and WorldBank datasets.

4. Adjusted Score for Invalid Data: this is the percentage of data that is not invalid for Mandatory Features (i.e. features without which the record would have to be discarded), simply speaking this is the count of cells containing valid data in mandatory columns Vs total number of cells expressed as a percentage. The higher this score the cleaner the dataset is from the perspective of invalid data. Validity checks are both generic and context specific as defined above. This score is available for both Starbucks and WorldBank datasets.

Validity checks implemented for Starbucks dataset

Validity checks implemented for Worldbank dataset

DQS for Starbucks and WorldBank datasets

The following code fragment downloads a CSV file containing these scores from the Github repo of this project. These scrores have been calculated offline by running the program and uploading the results as part of the Github repo. The missing data raw score for the WorldBank data is very less ~34%, this is because for several parameters the 2015 data is not yet available, so in this case we will use the 2014 data as these are macro level indicators which do not see a drastic change year on year (this will be disussed more in the next phase of the project).

Feature Creation

Three (for now) new features are created using the datasets. These features are all added to the Starbucks dataset. These are described below:

The following table shows a view of the modified dataset with new features added.

Some notes about the code and how to run the SBS program

This section provides information about how to run the code, how is the code organized amongst different modules and what are the output files generated.

How to run the code

When you unzip the zip file from Blackboard Github, you will see a folder called sb_study created which contains all the code, helper files and even an output folder (the contents of the output folder get rewritten on every run of the SBS program). As a pre-requisite the machine on which this program is being run is required to have Python 3.5.2|Anaconda 4.1.1 installed. Run the program by simply saying

python main.py

The code will take about 2 to 3 minutes to run completely. As the program runs it would spit out traces on the console and these traces are also automatically being logged into a file called SBS.log created in the output directory (also created by the SBS program in the current working directory if not already present). The program would create a directory called output and several sub-directories inside it to store the various artifacts created (csv files, plots etc.) as part of running the program.

The program can also be run in analysis only mode such that it assumes that the data has already been downloaded by running the program in the full mode first. In the analysis only mode the program runs various forms of analysis on the data and stores the created artificats under subdirectories in the output folder. This mode is useful if only the analysis needs to be redone and the entire datsset does not need to be downloaded again. To run the program in analysis mode use the following command line.

python main.py -a

Code organization and other files

Here is a listing of everything that is included alongwith a brief description.

A note about the structure of the output folder.

The output folder contains several csv file which correspond to data as downloaded, cleaned data and combined data.

• WDI_SB.csv: this is the combined dataset, all analysis is run on this dataset.
• DQS: this folder contains various dqs generated at different stages, contains multiple files, the contents are self-explnatory. The Jupyter notebooks displays the contents of these files.
• EDA: this folder ontains artifacts from the exploratory data analysis, which include csv files and plots, the contents are self-explnatory. The Jupyter notebooks displays the contents of these files.
• scatter: this folder contains various scatter plots generated for the combined dataset. Specifically, the plots between number of Starbucks store and the WDI indicators.
• regression: this folder contain various artifacts (csv files and plots) generated as part of regression.
• clustering: this folder contains various artifacts (csv files and plots) generated as part of clustering.
• classification: this folder is a placeholder for artificats created as part of classification. Currently the classification results can be seen in the Jupyter notebook and the SBS.log file.
• association_rules: this folder contains various artifacts (csv files only) generated as part of association rule mining. These csv files are included as part of the Jupyter notebook.

Combinining the two datasets

The WorldBank WDI indicators dataset and the Starbucks datasets are combined together for hypothesis testing and machine learning described later in this notebook. To combine the two datasets we first filter out all the countries from the WorldBank dataset which do not have a Starbucks store i.e. we keep only those countries in the combined dataset that have at least one Starbucks store. This combined dataset is then cleaned (described later in this document) and all analysis is run on it. There are 72 countries in the world where Starbucks is present so the combined dataset contains 72 rows and 48 columns (more on the column count later). The following fields computed from the Starbucks dataset are added to the combined dataset:

The shape of the combined dataset is (72,48). This includes 35 WDI (World Development Indicators) from the WorldBank datset.

For converting continous columns to categorical a quantile (percentile) based binning strategy is used.

Exploratory Data Analysis

This section describes various types of exploratory analysis run on the two datasets as well as on the combined dataset before getting into the hypothesis building and machine learning part.

Basic Statistical Analysis and data cleaning

This section provides basic statistical analysis for both the WorldBank dataset and the StarBucks datasets. The following boxplot provides a quick overview into the distribution of Starbucks stores across continents. Some quick observations from the box plot are described below:

1. There is a huge variation in the distribution of Starbucks stores in the continent of North America. This is on expected lines since other than the U.S. and Canada all other countries are islands. Also the North American continent contains the country with the most number of Starbucks stores (the U.S. ofcourse, no surprises there).

2. There is minimum variation in the distribution of Starbucks stores in the continent of Oceania since it consists of only two Countries which have Starbucks stores i.e. Australia and New Zealend.

3. The European continent has the most even distributed around the median.

Mean, mode, median, standard deviation

The WorldBank dataset consists of all numeric values except for the country code, country name and the derived feature for ease of doing business (categorical). The country code and name are excluded from the statistical analysis since there is onnly one row per country in the dataset. The statistical measures for each feature are listed below. Some of the features show NaN for the statistical measures, that is because there is no data at all for these features (they are subsequently removed from the dataset, described later).

The Starbucks dataset consists of all categorical features, therefore only a mode calculation is valid. Several features such as street names are excluded from the mode calculation as well. There are some intresting observations to be made

1. The maximum number of Starbucks store in a single city are in a city outside of the U.S.
2. The maximum number of Starbucks store in the U.S. and in the Eastern Timezone.
3. The maximum number of Starbucks stores (about 70%) occur in countries with very high ease of doing business.

Outlier Detection

Outlier detection is performed on the WorldBank dataset. We determine outliers using a simple rule, anything outside of the 3rd standard deviation is an outlier. The results of outlier detection are provided below. Note that while there are values which get categorized as outliers but they are still valid values (as understood by human inspection of the data) therefore no outlier handling as such (smoothning, removal etc) is done for these values. Also given that these are valid values they cannot be ignored as they may have a significant role in the machine learning algorithms.

No outlier detection is performed on the Starbucks dataset as it contains only categorical values (values outside of well defined valid values set are considered invalid and that analysis has already been performed as part of the data quality detection phase).

Bin the data

As part of the analysis that follows several continous/discrete features were binned (as described below this is needed for association rule mining, feature selection etc.). Here we show one specific feature IC.BUS.EASE.XQ which quantifies the ease of doing business in a country (1 being highest, larger the value lower the ease of doing business) into a categorical value. Numeric values get mapped to categories as follows:

The benefits of binning this feature is for exploratory data analysis to answer questions like what % of Starbucks store are in countries with high or very high ease of business. The method used provides a simple strategy to bin the data, although it would require a subject matter expert to say what the most appropriate bin edges are.

Here is an excerpt from the WorldBank data showing the new column created as a result of binning along with the original feature (first 5 rows shown).

Handling missing values

There are not a lot of missing values in the Starbucks dataset (the adjusted data quality score for missing values is a 100% which means all the featurs that we care about have no missing values, see section on DQS above). For the WorldBank dataset has a lot of missing though, DQS for missing values is ~ 33.65%. The missing values in the WorldBank dataset are handled as follows:

1. The WorldBank data is updated every year and many times there is a timelag in the data being available for a year. Also since these indicators are macro level indicators so the year on year change is not drastic. As a first step, all the values that are missing for the latest year are filled with the values from the previous year whenever the previous year's data is applicable.

2. Define a feature density threshold of 90% and delete all features which are less than 90% filled.

3. Finally the analysis is to be performed on the combined dataset, as described above the combined dataset is obtained by joining the WorldBank and Starbucks datasets only for those countries that exist in the Starbucks dataset. This intersection results in a dataset which is pretty full (99.8%).

   • Any missing values are handled by replacing them with the mean value of the feature.

The DQS metrics for the original datasets as well as the combined dataset after cleaning are shown below.

Data cleaning for the combined dataset

The combined dataset is cleaned for missing values as described in the previous section.

1. Missing values are handled as described above (take missing data from 2014, remove features for which density is less than 90%, many rows get removed when we combined the WorldBank data with Starbucks data).
2. There are no invalid values in the WorldBank dataset (see Data Issues and Data Quality Score sections above). The invalid values from the Starbucks dataset (some lat/long values were missing, names of some streets could not be read due to errors with unicode handling) do not impact the overall analysis and are therefore ignored, see adjusted DQS for invalid data.

Histograms and Correlations

This section provides histograms and corelations from the individual datasets where appropriate and also from the combined dataset used for analysis.

Histograms from WorldBank dataset

The following figure shows histogram of three different features from the WorldBank dataset.

The observations from these histograms are described below.

Histograms from the combined dataset

The following figure shows histogram of three different features from the Combined dataset. The value of the coorelation coefficent is also shown.

The observations from these histograms are described below.

Scatter plots and correlation coefficient WorldBank dataset

The following figure shows the scatter plot for three features in the WorldBank dataset.

From the scatter plots we can observe that there is a negative correlation between ease of doing business and number of Internet users in a country, but, since the ease of doing business is expressed in a way that lower the number greater the ease so therefore it follows that ease of doing business and number of Internet users are positively correlated.

Ease doing business is also positvely correlated to per capita GDP. There also seems to be a positive correlation between number of internet users and per capita GDP.

It is important to mention here, the oft repeated phrase correlation does not imply causation.

The following table shows the correlation coefficient between these features and the value of r corrobrates the observations from the scatter plot above.

Scatter plots and correlation coefficient the combined dataset

The following figure shows the scatter plot for three features in the WorldBank dataset.

From the scatter plots we can observe Number of Starbucks stores has somewhat of a positive non-linear relationship with the number of international tourist arrivals, total population and number of people having access to the Internet. This is further clarified when we do polynomial regression between number of Starbucks stores and number of international tourist arrivals.

The following table shows the correlation coefficient between these features and the value of r corrobrates the observations from the scatter plot above.

Cluster Analysis

Cluster analysis was conducted to find patterns within data (from the combined dataset). Recall that the combined dataset is created by filtering out countries which do not have Starbucks presence, this brings a Starbucks specific dimension to this dataset. The presence of Starbucks stores (lets say in terms of number of stores in a country or any other feature derived from the Starbucks dataset) is the dependant variable in this analysis, therefore any patterns we find in the clusters that are created would be examined in that light i.e. do the cluster member show an affinity to a particular value of the dependant variable.

Before doing clustering a Principal Component Analysis (PCA) is first conducted to transform the data into 3 dimensions. PCA in 3 dimensions helps because the clusters can the be visualized in 3D space. Three different types of clustering mechanisms are used viz. KMeans, DBSCAN and Hierarchical (Agglomerative) clustering. This is described in detail in the next few sections, lets examine the PCA results first.

The following figure shows the result of PCA in 3 dimensions. Note that PCA only helps to spatially distinguish the data, but seeing the clusters without a spectral differentiation is difficult (which is what is provided by the clustering algorthms).

KMeans, DBSCAN and Heirarchical Clustering

The following figures show the result of 3 clustering techniques. It appeas that KMeans(k=5) provides the best custering.

Observations from clustering

1. The plots above indicate that KMeans has the best clusters, although it has to be said that the clusters are not very clear and spaced out even with KMeans which indicates that there may not be very clear patterns in the data.
2. Value of k=5 was choosen as it matches the bins that we want to create for classifying number of Starbucks stores as Very Low, Low, Medium, High and Very high.
3. Vectors at the center of the clusters (KMeans) are centroids and the points in a cluster are all placed around it.
4. The hierarchical and DBSCAN clustering produce similar results as can be observed just by visual inspection of the plots, the KMeans clusters seem to be much more well formed.

Specific observation from KMeans

The following table shows a filtered list of labels from each of the clustering algorithms attached to each entry in the dataset. The KMeans labels does seem to show a pattern in countries that are grouped together.

1. KMeans label=4 groups together U.S. and China in a cluster of their own,which should not be a surprise. Ofcourse we can see any number of similarities that could have caused this for example GDP, tourist arrival etc. Also from a Starbucks perspective, these countries represent very high number of Starbucks stores (top 2).

2. KMeans label=3 groups together Brazil, Russia and other countries of the erstwhile Russian federation. From a Starbucks perspective these countries represent very similar Starbucks store density (number of stores per 100,000 people).

3. KMeans label=2 groups together mostly Asian, African and some South American countries which have low to medium ease of business but medium to high population.

Other labels can also be explained in a similar way.

Association Rules / Frequent Itemset Mining Analysis

Association rule mining is run on the combined dataset to find tuples of various feature values that occur together in the dataset and draw inferences from that. Most of the data in the combined dataset is continous or discrete, to run association rule mining the first step that is needed is to convert data into categorical values. The combined dataset contains 40 odd features, so instead of converting the entire dataset to categorical, some selected features are converted and association rule mining is run for those features.

Continous/discrete features are converted into categorical by binning them by percentile values as described in previous sections. The dataset for this particular problem does not lend itself very well for this type of analysis (because large number of features are numeric). The conversion of numeric data to categorical in the context of this dataset is an art rather than a science (meaning what bin sizes are good for categorizing number of stores as low, medium or high, this is subjective) therefore we can get vastly different results by changing the criteria for converting numeric data to categorical (have 3 bins instead of 5 for example). Regression and classification algorithms discussed later are better suited for this data science problem.

To select which features would have a relationship with dependant variables that we want to predict we draw scatter plots of the dependant feature Vs all numeric features in the dataset and manually examine where a correlation is seen. Another technique used for finding out the most important feature using a machine learning technique (ExtraForestClassifier) and then using the top two most important features for conversion to categorical. The scatter plots are not being included in this report allthough they are available as part of the overall package that has been submitted. Similarly the list of important features can be seen as part of the logs generated on running the program.

Two different types of rule mining is done, one with a rule containing a single variable predicting the dependant variable and another one with two variables predicting the dependant variable.

Association rules with 2 features

The following table identifies all th association rules with two variables:

The expectation is that very high influx of tourists should occur together with very high number of Starbucks stores. This does occur as this is the second most frequent rule as shown in the table below, however the support value for this rule is a low 0.125 and the confidence is 0.64. The rule that occurs most frequently is medium tourist influx corresponds to medium Starbucks stores, this rule occurs 15 times and has a support of 0.20 and confidence of 0.53. Expected to have a greater support level for the more frequent rules but this not happening probably because of two reasons a) there are only 72 countries with Starbucks stores so there is not a lot data and b) the conversion of numeric to cateogrical data in this case is not an exact science and can be tweaked to give better results (reduce the number of categories for example).

This is all summarized in the following table. The rules themselves are written in the following format R:({A})->B means if A is present then B is also present (with the frequency, support and confidence mentioned alongside to give a sense of probability).

Association rules with 3 features

The following table identifies all th association rules with two variables:

The expectation is that high influx of tourists and a high local population should occur together with high number of Starbucks stores. This does occur as this is the second most frequent rule as shown in the table below, however the support value for this rule is a very low 0.08 but the confidence is 0.75. The rule that occurs most frequently is medium tourist influx in a country with medium population corresponds to medium Starbucks stores, this rule occurs 7 times and has a support of 0.09 and confidence of 0.53. Expected to have a greater support level for the more frequent rules but this not happening probably because of two reasons a) there are only 72 countries with Starbucks stores so there is not a lot data and b) the conversion of numeric to cateogrical data in this case is not an exact science and can be tweaked to give better results (reduce the number of categories for example).

This is all summarized in the following table. The rules themselves are written in the following format R:({A,B})->C means if A and B present together does imply C (with the frequency, support and confidence mentioned alongside to give a sense of probability).

Association rules with 3 features filtered by 3 different support levels

The following tables present association rules filterd by 3 different support levels.

1. support >= 0.05 and confidence >= 0.25
2. support >= 0.07 and confidence >= 0.25
3. support >= 0.1 and confidence >= 0.25

Predictive Analysis

This section describes three hypothesis that have been developed using the combined WorldBank and Starbucks dataset and then tests each of the three hypothesis using parameteric and/or predictive algorithms.

Hypothesis 1: No difference in the average number of Starbucks stores in countries with very high number of international tourists Vs rest of the world.

This hypothesis is tested in two ways.

1. Two sample T-test: A two sample T-test is used to find out if the means of two independant samples are different. The null hypothesis is that means of both the samples are the same. The two sample t-test only relies on sample means and does not need to know the population mean.

2. Linear and Polynomial regression: ordinary least squares regression. Linear regression tries to predict the value of the dependant variable by using a linear combination of explanatory variables. Linear regression could be univariate (one explanatory variable) or multi variate. A polynomial regression tries to model the dependant variable using higher powers of the explanatory variable(s).

Using the two sample T-test for hypothesis 1

The combined dataset is used to filter out entries belonging to 'ST.INT.ARVL.Categorical' equal to 'VH' (for very high) and then another set of entries corresponding to 'ST.INT.ARVL.Categorical' value other than 'VH' so as to include the rest of the world. Python scipy module is used to calculate a T statistic and a corresponding p-value when comparing the two distributions. The p-value comes out to be ~0.17. This is a unusually large value and with a typical T-test alpha value of 0.05 it would not be possible to reject the null hypothesis but in this case since the distribution of the number of stores was assumed to be normal when it is known that it is not normal (see histogram in previous sections) so a higher alpha of 0.20 is choosen and this enables the null hypothesis that there is no difference to be rejected. This is in line with what is seen in the scatter plots for these two features.

Using Linear and Polynomial regression for hypothesis 1

In this method instead of using the categorical value of the number of international tourist arrivals the actual discrete value is used and it provided to a linear regression model and a polynmial regression model to predict the value of the number of Starbucks stores. The linear model provides an explained variance of 0.27 whereas the polynomial model performs much better by providing an exaplined variance of 0.33. Ideally the the explained variance should be close to 1. It is a matter of further investigation as to what else can be done to make the prediction better. The following scatter plots show linear and polynomial regression models in action, as can be seen from the plots the polynomial model provides a better fit. The polynomial model is a degree 3 polynomial (other degrees were also tried and 3 provided the best results). The outlier value seen towards the stop right is the number of stores in te U.S., it is an extreme value and further investigation would need to be done to figure out if a single variable can predict such a huge jump in value.

Hypothesis 2: The number of Startbucks stores in a country as a categorical value can be predicted using some WDI indicators

This hypothesis states that the number of Starbucks stores in a country expressed as a categorical value (Very High, High, Medium, Low, Very Low) can be predicted using WDI (World Development Indictors) from the WorldBank data. This is verified by running various classification algorithms and then selecting the on that provides the highest accuracy value as measured via cross validation.

The combined dataset contains 35 WDI features. All of these are not used for classification, the most important 15 features are choosen. The Extra Trees classifier which uses an ensemble of decision trees is used to provide a ranking for all the features and then the top ranked 15 features are used as inputs for the classification algorithms. Tests were done with 5, 10, 15, 30 and all 35 features as well, 15 most important features provided the best results in terms of accuracy.

The 15 most important features for determining the number of stores value as a categorical are:

The above list is a very intresting list because all of these indicators seem to appear as very important factors that a business would consider. A lot of these are either economic or directly related to business activity so it intutively makes sense that the machine learning algorithm was choosing the right indicators.

Once the indicators were selected all of the following classification algorithms were applied on the data and the results in terms of accuracy and standard deviation of a K-fold cross validation was noted. Random forest turned out to be the best classifier in terms of accuracy.

The accuracy turned out to be in the 50 to 60 percent range but this is just the first pass at predicting. The indicators, how they are used, and finally the bins for the classification of the dependant variable, all are subject to be refined and this is not the final score. The hypothesis is considered valid although the accuracy value is lower than desirable. More investigation needs to be performed to bring the accuracy score up around the 80 percent range.

Description of classification algorithms used

All of the following algorithms were used for classification. All of these algorithms are available via the Scikit learn python package. A test train split was done using training data fraction as 20%. A 10-fold cross validation was done to determine the accuracy of various algorithms.

Hypothesis 3: The number of Starbucks stores (discrete) in a country can be predicted using some WDI indicators.

This hypothesis stats that number of Starbucks stores in a country which is a discrete variable (as opposed to the categorical considered in hypothesis 2) can be predicted using some WDI indicators. The method used here is similar to hypothesis 2 in that first a set of three most important features are determined using the Extra Trees classifier and then these three features are fed into a Multivariate Polynomial Regression model to predict the dependant variable i.e. number of Starbucks store in a country. The three features used for this prediction are 'ST.INT.ARVL' (international tourist arrivals), 'TX.VAL.TECH.CD' (High-technology exports (current US dollars)), 'SP.POP.TOTL' (population of the country).

A multivariate polynomial model is choosen since there are more than one explanatory variables and the relationship with the dependant variable is definitely not linear as is known from the scatter plots (not included here but available as part of the package). A degree 3 polynomial works best as determined by experimenting with other values as well. An explained variance score of 0.99 is achieved, however the mean squared error is still very high at 27063. The MSE needs to be improved before the hypothesis can be considered valid, this needs further investigation in terms of feature selection and model used for prediction.

Other Hypothesis

A couple of other hypothesis were also tested as part of this phase.

1. Whether ownership type of Starbucks store ina country should be mixed (Starbucks owned, franchisee owned or other) or just Starbucks owned. All classification models were tried and it was found that Random Forest and kNN both predict this with about 84% accuracy (stddev 12%).

Bringing it all together [Project1 + Project2]

At this time in the data science lifecycle a clean combined dataset exists and has been used to validate multiple hypothesis as described above. The focus of this phase was getting everything in place to run predictive algorithms, the next phase will focus on improving the accuracy of the predictions. The results of the predictions are included in the overall package (see regression folder) but are not being disucssed here at this time.

What does appear is that more than a regression problem the problem can be modeled as more of a clustering and classification problem. In other words, use clustering algorithms (like KMeans as described above) to see if countries which have no Starbucks stores appear in which clusters, for example if they appear in clusters corresponding say countries with a high number of Starbucks store then that is a prediction. Similarly, the RandomForest or NaiveBayes classifier (which provided between 50 to 61% accuracy) could be used to predict the category (low, medium, high etc) for number of Starbucks stores.

Another aspect to review is that maybe the categories for the number of Starbucks stores need to be reduced from 5 (Very high, high, medium, low, very low) to just 3 (high, medium, low) to improve the accuracy of predictions. These are all things to be considered for the next phase of the project.

Predictions based on Classification

As discussed in the previous section, it appears that classification is a better option for prediction than regression for this data science problem. In other words, it is possible to predict the number of Starbucks stores in a country as a categorical variable (i.e. VeryLow, Low, Medium, High, Very High) with a much higher degree of accuracy than predicting the exact number as a continous variable. The RandomForests classifier discussed above provided an accuracy of between 50 to 61% (stddev 12 to 20%), it is used to predict the number of stores in countries which do not have a Starbucks presence yet but for which there is data available for all the features to enable prediction. The features used for prediction are listed in Hypothesis 2 and are not being repeated here. Since the original WorldBank data has a lot of missing data therefore predictions cannot be made for all countries where Starbucks does not have a presence yet. The following table lists the results of the prediction. (VL: Very Low, L: Low, M: Medium, H: High, VH; Very High). The categories translate into the following numerical values:

Discussion of the prediction results

The prediction throws up some interesting results. While some interpretation is being provided here, it would upto a subject matter expert (and well, Starbucks) to say why these groupings make sense (or not). A point to note here (which has not been discussed earlier) is that the WolrdBank data also had data for groups of countries (such as OECD) combined as a single entity so these groups also appear in the following classification.

Time Series Analysis for number of Starbucks stores [Additional Analysis 1]

A piece of information available in the Starbucks dataset is the date on which a store was first seen and the dataset contains information from December 2013 to the current month i.e. November 2016. The first seen field can be used to identify (say) a month on month increase in the number of stores i.e. stores for which the first seen date is say a date in January 2014 can all be considered as newly opened stores in January 2014 (this was also empirically verified by the fact that a new store that opened in Clarksburg, Maryland that the author is personally aware of shows up with a first seen date of November 2016).

This information about the first seen date can be used to do a time series analysis to examine and predict the number of Starbucks stores in the next 12 months. To do this the following steps were followed:

1. Extract a time series of interest i.e. total number of Starbucks stores across the world and plot it to visually examine how the number of stores was increasing over time.

2. Examine the time series for being stationary i.e. is the mean of the series increasing with time. The idea behind requiring the time series to be stationary is that if the time series is stationary then individual points can be considered I.I.D (independent identically distributed) and the usual concepts of stationary random variables (Central Limit Theorem, Law of Large Numbers) can be applied to it. If the series is not stationary then it needs to be made (read transformed into a) stationary series. There are several transformations to stationarize a series such as Logarithmic, Deflation by CPI, first difference etc.

   • From visual examination of the number of Starbucks stores time series it is clear that the series is not stationary (number of stores is increasing with time). A Dickey-Fuller test could also have been done to check for stationarity but it has not been done as part of this project. The logarithmic transformation was used to stationarize the time series.
   • Seasonal decomposition of the Starbucks stores time series is also done to extract out trend and seasonal variations. Plots are presented below.

3. The Auto Regressive Integrated Moving Average (ARIMA) model is used to forecast the timeseries. In general, ARIMA models are applied in some cases where data show evidence of non-stationarity (like in case of the Starbucks store data), where an initial differencing step (corresponding to the "integrated" part of the model) can be applied to reduce the non-stationarity. See https://www.otexts.org/fpp/8/1 for more details about time series and https://people.duke.edu/~rnau/411arim.htm for an explanation on ARIMA.

   • The parameters (p,d,q) of the non-seasonal ARIMA model are are non-negative integers, p is the order (number of time lags) of the autoregressive model, d is the degree of differencing (the number of times the data have had past values subtracted), and q is the order of the moving-average model. The (p,d,q) values are generally identified by studying the ACF/PACF values, in this case, they were identified by trial and error (the values that gave the least mean square error were used).

   • A Seasonal ARIMA (SARIMA) could have been tried but current stable release of the statsmodels package (0.6.1) does not support it, although it is available in a version currently under development (0.7.1).

4. The ARIMA model is provided the logged timeseries for reasons discussed above. The fitted values provided by the model are actually the differences of the logged series, these values are then cumulatively summed to derive the actual logged series. Once the predicted logged series is obtained then it is transformed back to its original value by taking an exponent (e^x) of the logged value. The forecasted values for future timestamps is obtained using the same procedure.

The following sections discuss time series analysis for two timeseries, one for the total number of stores across the world and another one for the number of stores in the state of New York. Several other TSA analysis (such as time series for number of stores in the U.K, India, the continent of Africa) are included in the project submission.

TSA for Starbucks Stores across the World

There are four plots provided below.

1. Seasonal decomposition: clearly shows that there is an almost linear upwards trend as well as a seasonal element (a frequency value of 5 was used).
2. First difference of the logged difference along with the fitted value from the model.
3. Original series with the fitted values (transformed back from the logged differences provided by the model). Also shown here is the root mean squared error (RMSE) to determine how good the prediction is. It is seen that the RMSE is ~152 which is roughly 0.60% of the total number of stores (25090), so the fitted value is off by less than 1% of the original value.
4. Original series with fitted value and projected values. The projected values are for the next 12 month period i.e. from January 2017 to December 2017. The projected number of Starbucks stores in the World at the end of 2017 should be around 27995.

TSA for Starbucks Stores in the state of New York

Same four plots as above are provided for the timer series analysis for number of stores in the state of New York as well. The RMSE is again less than 1%. Forecasted number of stores in the state of New York at the end of 2017 is ~692, which would be an addition of 53 stores from the current number of 639.

Exploring distribution of Starbucks stores across the world [Additional Analysis 2]

This analysis explores the distribution of the number of Starbucks stores in countries which have a Very High density of Starbucks stores (> 1 Starbucks store per 100,000 people) and at least 100 Starbucks stores across the country. The number of stores per city in the country is modeled as a random variable and its distribution is examined. The following table provides the list of countries that have been studied, note that Singapore which has a 125 stores across the country is not included in this list since there are only two Cities in Singapore one of which has a single Starbucks store.

From the above table it is clear that the maximum number of stores occurs in the captial city of the country except in the case of the U.S. where the maximum number of stores in a single city is in New York (which is not a surprise).

The following boxplots and histogram show the distribution of Starbucks stores across cities in the above mentioned countries. As is clear from these plots, the distribution is right skewed with very large difference between the third quartile and the maximum (because the one city with the maximum number of Stores is an outlier).

In order to obtain any insight from the above plots the outliers need to be removed. This is done by "Winsorizing" (see https://en.wikipedia.org/wiki/Winsorizing) the distribution at 0.05 level i.e. all greater than the 95th percentile are clamped at the 95th percentile value. The following are the plots after Winsorization. As seen from the plots it is evident that the basic nature of the plots has not changed after Winsorization but at least now it is clear that median number of stores per city for all of these countries is around 3 and the 3rd quartile is less than 7. The shape of the distribution as seen from the histograms is similar for all of these countries with the U.S. having the maximum number of stores in almost all percentiles. The shape of the distribution is discussed in the next section.

Looking at the histogram above, this distribution looks similar to an exponential distribution or at least something from that family of distributions. To understand this better, the distributions for each of these countries was examined individually and it was found by visual inspection that a Exponential Logarithmic distribution or an Erlang distribution fitted this best. The Seaborne package for making plots also provides the option to define a "fit" function and then it overlays the fitted plot on the actual data.

The following plots show the CDF (cumulative distribution function) plots as well as the histograms with the fitted curve from the Erlang distribution. The plots provided here are for the U.S. and the U.K. (rest of the plots are included in the submission package). The Erlang distribution is used to model incoming calls in a telephone system, waiting times in queuing , business economics for describing interpurchase times etc. In the context of this data science problem, we could think of applying this distribution to model number of Starbucks stores in cities of a country where Starbucks does not exist today. The distribution could define the number of stores in the capital city (which would be an outlier) and number of stores in other cities in the country which would be within the 1st quartile, 3rd quartile and so on. This could be used in conjunction with the previous analysis where the number of Starbucks stores is predicted as a categorical variable, so from there a range for the total number of stores in the country is known and this range could then be modeled using the Erlang distribution to distribute the total number of stores across multiple cities. This is a topic for further research.

Visualizations

This sections provides several visualizations done as part of this project to better understand the data at hand. These visualizations are discussed individually in the subsequent sections.

Starbucks around the world [Visualization 1, Bokeh]

The following visualization shows all Starbucks stores as of November 2016, overlayed over a world map. The brown dots on the map represent locations of Starbucks stores. A hover over individual point brings a popout text which provides the name of location where the store exists along with the name of the city. Since the map is displayed from a certain zoom level therefore the points representing Starbucks stores appear very close to each other in high density areas.

The visualization clearly shows that a lot of Starbucks stores are located on the eastern and western coast of the US, the UK in Europe and China and Japan in Asia as well as Mayasia and Indonesia in south east Asia. There are relatively fewer stores in the South America and extremely few stores in Africa.

This visualization has been created using GMapPlot and Bokeh.

Starbucks in the U.S. [Visualization 2, SVG]

With more than 13,000 Starbucks stores, the U.S. has the most Starbucks stores in the World by a huge margin, but how are Starbucks stores distributed within the US? This is what this next visualization tries to explore by creating a heatmap of U.S. counties, more the number of stores, more intense the color of the map. Doing a mouse over a county brings up the name of the county and the number of stores in that county. The stores are expectedly most on the east coast and the west coast and very few in the states at the geographical center of the country like in Nebraska, South Dakota, Kansas etc.

Maximum number of stores, 665, occur in Los Angeles county in California. New York county (234) in New York state, Maricopa county (334 stores) in Arizona, Cook county (325 stores) in Illinois, King county (336 stores) in Washington are some of the counties with a very high number of Starbucks stores. Just from visual examination of the graphic it appears that the Ney York county (which includes Manhattan) would probably have the highest density of Starbucks stores in terms of number of stores divided by the size of the county.

This plot has been created using a combination of an SVG (Support Vector Graphics) image of the US map with counties marked on it and the FCC (Federal Communications Commission) API that provides the name of the county given a locations latitude and longitude. All Starbucks locations are mapped to a county name using the FCC API in an offline process and a count is maintained for the number of stores in each county. Once the store to county mapping is complete the result is stored in a CSV file. This CSV file is then used to assign color and hover text to the counties as listed in the SVG file. The counties are colored using shades in a single color pallete with 0 (no Starbucks stores) in the county being the lightest to greater than 300 stores being the darkest.

Starbucks stores | Internet Users Vs International Tourists Arrivals [Visualization 3, Plotly]

As discussed in previous sections, the number of international tourist arrivals per year and internet users (per 100 people) are two important features which have been seen to influence the number of Starbucks stores in a country. These features were determined to be two of the more important features by the ExtraTreesClassifier and even intuitively it does make sense that countries with a large number of tourist arrivals and high Internet accessibility should have a higher number of Starbucks stores.

This relationship between number of Starbucks stores, Internet users and International tourists is depicted using a bubble chart. The bubble chart represents 4 dimensional data as described below.

1. The X axis represent the international tourist arrivals per year.
2. The Y axis represent the number of Internet users per 100 people.
3. Each country is represented by a bubble, the diameter of the bubble representing the number of Starbucks stores in the country. Higher the number of stores, larger the bubble size.
4. The color of the bubble represents the continent to which the country corresponding to the bubble belongs to.

The bubble chart presented below is created using Plotly, mouse over the bubbles to see information about the country, number of Starbucks stores etc. As with all Plotly charts, basic user interactions like selecting specific areas of the chart, saving as png etc are available. All the data required for the bubble chart is available as part of the combined World Band Starbucks dataset.

Some important observations from the above visualization are as follows:

1. Almost all the bubbles which are medium to large in size all occur towards the right side of the 10 Million mark on the x-axis which indicates that unless a country has north of 10 Million international tourists arrival a year it would have very small number of Starbucks store probably around 20 or so with the only exceptions being Brazil, Argentina, Philippines and Indonesia.

2. The green colored bubbles which correspond to Asia are very visible (most of them are at least medium size with China being the largest). This seems indicates that countries in Asian generally tend to have a higher number of Starbucks stores than counties in other continents. This is also indicated in the boxplot drawn in a previous section, the mean number of Starbucks stores was the highest for Asia.

3. The largest bubble towards the top right side of the chart is ofcourse the U.S., high availability of Internet access and a very large number of tourists.

4. France which is the blue colored bubble on the top right corner is relatively much smaller than bubbles representing other countries with high availability of Internet access and a large influx of tourists. This looks like an aberration, France has only 124 Starbucks stores, at more than 80 Million tourists per year and more than 80 out of 100 people having access to the Internet, France should be having a lot more Starbucks stores than just 124.

Increase in Starbucks stores worldwide [Visualizations 4 & 5, Plotly]

The following bar chart shows the increase in Starbucks stores by country from December 2013 to November 2016. Mouseover individual bars to see the name of the country and the number of Starbucks stores added in that country between December 2013 to November 2016. As with all Plotly charts, basic user interactions like selecting specific areas of the chart, saving as png etc are available.

Some observations from the above visualization:

1. The maximum increase in terms of absolute number of stores happened in the U.S. followed by China, South Korea, Canada and Great Britain.
2. Most of the growth is concentrated in the top 10 countries, in fact the increase in the number of stores in U.S. and China put together is way more than the increase in the number of stores in the rest of the world put together (almost 1.35 times more).

While the increase in absolute terms is useful, it is also interesting to find out which countries have the highest rate of increase, basically where is Starbucks expanding fastest. This list may not be the same as the list of countries which have the greatest growth in terms of absolute numbers. To determine this, the following procedure is follows:

1. A month on month increase in number of Starbucks store is determined for each country.
2. It is then converted to a percentage of the next month's value.
3. This percentage month on month increase is then totaled up for the entire duration from December 2013 to November 2016.
4. The countries with the highest cumulative percentage increase are determined.
   • Note that countries with less than 25 stores as of November 2016 are excluded from this analysis (anything less than 25 stores is considered too small to be significant).

The timeseries for the number of stores for 15 countries having the highest cumulative percentage increase is represented below as Sparkline charts (more on Sparklines https://en.wikipedia.org/wiki/Sparkline). Note that this list does not include the U.S. so even though the U.S. had the highest increase in terms of absolute numbers the fastest rate of increase is somewhere else. It is seen from this visualization (and confirmed from the actual numbers, not shown here) that the maximum cumulative percentage increase happens in India where the number of stores increased from 29 in December 2013 to 85 in November 2016.

Limitations and future work

The following limitations and scope for future study is identified:

1. Many WDI features had to be discarded because there just wasn't sufficient data available, it is possible that some of these missing features might have been important from a prediction perspective. As and when data for these features becomes available then those features can also be included in the analysis.

• Predictions could not be provided for all countries where Starbucks does not have a presence because data for features used for prediction was not available for these countries, once this data becomes available then these predictions can also be provided.

1. Are we able to increase the accuracy of predictions by reducing the number of categories to 3 (Low, Medium, High) instead of the current 5? This will create a tradeoff between having these categories defined over much narrower interval (and hence being more effective) and the accuracy with which these categories can be predicted.
2. The python statsmodels module would provide SARIMA (Seasonal ARIMA) at some point (hopefully soon) then the TSA analysis would be done with that model as well.

• The parameters of the SARIMA module would be determined by observing the ACF/PACF values, something which has not been done in the current analysis.

1. The distribution of Starbucks stores across geographies should also be examined by taking a log of the count of Starbucks stores in a city and then creating a histogram would transform the distribution into something more manageable. See plot below. As with all Plotly charts, basic user interactions like selecting specific areas of the chart, saving as png etc are available.