## [R][Visualization] Radar Plot with Scotch whisky data

When  I had an Islay single malt for the first time, it was mind blowing. In my first foray into the world of whiskies, I took the plunge into the smokiest, peatiest beast of them all — Laphroig . The smell wasn’t pleasant initially but it got totally different when I took a sip from the glass. That same night, dreams of owning a smoker were replaced by the desire to roam the landscape of smoky single malts. What even dragging me more is how the same whisky can taste different as it ages.

As a relatively new scotch whisky fan, I wanted to investigate whether distilleries within a given region do in fact share taste characteristics. For this, I used a dataset profiling 86 distilleries based on 12 flavor categories.

From math&Statistics department at University of Strathclyde, a professor who seems passionate about whisky created a dataset with 86 scotch whisky distillery. This data includes 12 different taste categories, postcode and latitude and longitude. In this post, I will focus on creating radar plot for each whisky to show the range of tastes.

Data Description

Previewing the data,

86 malt whiskies are scored between 0-4 for 12 different taste categories including sweetness, smoky, nutty etc.

Then, I subsetted the data excluding unnecessary information for this post.

Required Library

Here are libraries that we need. For this post, ggRadar function is the main one and it’s from ggiraphExtra library.

Main Code

For this case, I just selected the whiskys I know. You can use sample function in R if you want to see random whisky taste plots.  As you can see, ggRadar function is pretty straight-forward. Since we want to see taste profile for each distillery, let’s set ‘group =distillery’.

Tara,

Now with this code, we, the passionate whisky explorers, will easily identify the flavors of whisky and explore different kinds!

## [R][Visualization] Animation in R

Recently, I found there is a sub-reddit for R Visualization and ran into animation visualization. Through the reddit, I happened to check this blog, revolutionary analytics showing the animation of predicted Japan’s aging population from 2010 to 2050. When I examined each line in the blog, the package named ‘animation’ suddenly attracted me and I did more research on it for a couple of days.

When I was in college, as a mathematical visualization nerd, my hobby was drawing arts in my fancy TI 89 calcualtor(2010-2018, RIP).  The very first artwork by calculator was this heart.

The amazing package ‘animation’ was created by Yihui Xie and his package is also simple to use. You can find his complex examples in his website by simply googling his name.

In order to create the good looking heart curve, I refer Wolfram Alpha formula.  This code will generate GIF image of heart but you can also create HTML and PDF. In this code, I first created X and Y matrices which indicate the space of x axis and y axis in this visualization. By applying ani.option, it controls the behavior of the animation and this example restricts  the maximum number of steps in a loop to 10. Then, the function created the loop and the heart curve formula is nested in the for loop. In the for loop, you can see that a is the sequence from 0.9 and 1.3 and it makes the heart shape shrinking.

Then, here is the final product. Enjoy!

## [R][Statistics]Time Series Analysis with S&P 500 Stock price

I have not been writing statistical/ ML blog posts for a while. So, it’s good to come back!

Today, I will demonstrate how to apply time series analysis on forecasting stock market price. I won’t go over deep theory of time series analysis but will show the most fundamental model of time series analysis model.

1. Open the libraries that we need. Then, open the data.

2. Exploring the data

Here are the first 6 rows of the data. It contains date, open, high, low,close price and volume. The last variable indicates the initial of stock market. In this example, AAL is American Airlines Group.

By using summary command, it’s easy way to get an overview of the data. In this case, there are a few NA values in price variables. Getting rid of NA values would make data handling easier.

The prices approximately ranges from 1.5 to 2067 while the median values are around 62-63. It indicates that the price distribution is likely to be skewed. We can check it visually soon.

By using str command, we can check each variable’s format. We can see that date variable is treated as a factor variable. It would be better to change the data variable into date format.

3. Data Cleaning

From the previous step, 1. there are NA values that we want to discard, 2. transforming the date variable format.

Just getting rid of them would cause the difference in each variable length.The good news is that there were not many NA values in each price variable. Simply, we can change NA value to 0.

Then, we had trouble with the format of data variable. We can change it into the right format by using as.Date command.

Let’s double check how the data is altered.

Now, there are no NA values. Instead, the minimum values are all 0 across open, high,low and volume. Date variable is now in proper date format.

4. Histogram-Check the distribution of price variables

The price distributions are quite similar to each other. As we expected, the distributions are kind of right skewed which means the median is lower than the mean values.

5. Time Series Analysis Model- A bit of Theory

The most basic time series models are AR,MA and ARMA model.

AR(Auto Regressive Model)

Autoregressive (AR) models are models where the value of variable in one period is related to the values in the previous period. In other words, Autoregression is a time series model that uses observations from previous time steps as input to a regression equation to predict the value at the next time step.  AR(p) is a Autoregressive model with p lags.

MA(Moving Average Model)

Moving average (MA) model accounts for the possibility of a relationship between a variable and the residual from the previous period. MA(q) is a Moving Average model with q lags.

The role of the random shocks in the MA model differs from their role in the autoregressive (AR) model in two ways. First, they are propagated to future values of the time series directly. Second, in the MA model, a shock affects X values only for the current period and q periods into the future; In contrast, in the AR model, a shock affects X values infinitely far into the future

ARMA

ARMA is the combined version between AR and MA.  The AR part involves regressing the variable on its own lagged (i.e., past) values. The MA part involves modeling the error term as a linear combination of error terms occurring contemporaneously and at various times in the past. It’s denoted as ARMA(p,q).

Assumption for these models:  The variance is constant while the mean fluctuates

6. Create the function that creates time series object and Plot it

tsclean() is a convenient method for outlier removal and inputing missing values

ts() is used to create time-series objects

In this example, I used 3M stock(MMM) . Recently, the stock price plunged right after the most recent earning call.
7. Stationarity- A little bit of more theory
A stationary process has a mean and variance that do not change overtime and the process does not have trend. It’s a common assumption in many time series techniques is that the data are stationary.
But this time series does not look stationary. To confirm that we will use “Dickey-Fuller test” to determine stationarity.
In general sense,when the p value is smaller than 0.05, we can reject the null hypothesis. Since the p value is much higher than 0.05, we can’t reject the null hypothesis. In plain English, the series not stationary.
8.Differencing in time series
Differencing is a common solution used to make the variable stationary. Think about the concept of differentiation in basic calculus class.
After applying differencing, it looks much more stationary.
9. Select Candidate Model
In the step 5,  I mentioned the three time series model: AR(p),MA(q) and ARMA(p,q). Now this step will discuss which model to select and which lag(p or q) to pick. The most popular way to decide which model to use is examining ACF and PACF plot which are autocorrelation function and partial autocorrelation function.

The blue line above shows significantly different values than zero. Clearly, the first graph above has a cut off on PACF curve after the 2nd lag or the 3rd lag which means this is mostly an AR(2)/AR(3) process.  The second graph above has a cut off on ACF curve after the 1st or 2nd lag which means it will be MA(1) or MA(2). But at the right side of the graph, there is a lag above blue line. So, let’s try to use another method of selecting candidate.

R has function that automatically chooses the most suitable ARIMA model.

According to the result, ARIMA(2,1,2) would work the best.

10. Forecasting

At least, the forecasting said the stock market price of 3M is likely to be increasing. So, good news for 3M investors!

Potential Next part

Using Recurrent Neural Network

## [R]Network Analysis with Star Wars

Some of you(including me) must have been excited about the upcoming Star Wars Sequel movie this December. It’s less than a couple of months away until the released date! To celebrate it, I worked with Star Wars data network analysis.

I found the data from github and the data only contains characters in “Star Wars 4: A New Hope“.

1. Let’s call the data we need

The first step is to read the list of edges and nodes in this network. Those are fundamental elements for network analysis.  What are edges and nodes?  Let me explain this way. The World Wide Web is a huge network where the pages are nodes and links are the edges.  Visually, nodes form circles while edges form directions in network analysis.

In edges data, for example, the first row means there were 17 scenes when C-3PO and R2-D2 were together.

2. Call the library we need and form the data frame for network analysis by assigning edges and nodes.

What does it mean? – `U` means undirected
– `N` means named graph
– `W` means weighted graph
– `22` is the number of nodes
– `60` is the number of edges
– `name (v/c)` means name is a node attribute and it’s a character
– `weight (e/n)` means weight is an edge attribute and it’s numeric

3. Plot the graph

This is the simplest way of drawing network plot. However, it doesn’t look neat. Let’s add more options.

We can see that R2-D2,C-3PO,Luke,Leia, Han, Chewbacca and Obiwan are at the center. In other words, they are the center characters in Star wars 4:A New Hope. ALso, we can clearly notice that Darth Vader, Motti and Tarkin are forming a group and it indicates that they are likely in the same group(dark side).

But if you want to know the importance of the characters, this graph does not explain enough. Let’s do some extra works.

Now we cans see that characters located at the center have bigger circles that characters located in peripheral areas in the graph. `strength` will correspond to the number of scenes they appear in. And we’re only going to show the labels of character that appear in 10 or more scenes.

It would be interesting to see the peer groups among characters.

Most networks have a single giant connected component that includes most nodes. Most studies of networks actually focus on the giant component. `igraph` also makes it very easy to plot the resulting communities.

It is the simple version of network analysis. Hope you enjoy reading this!

Hello, for this post, I will show how to visualize spatial data on Google Map using R. It is simpler than you think.

What is Spatial Data?

it is the data or information that identifies the geographic location of features and boundaries. The data that I’m using today has longitude variable and latitude variable so that we can locate the data points accurately on the map.

Now you know what spatial data is roughly so let’s jump into the map visualization.

In ggmap, you need ggplot2 package.

ggmap library contains all the information of google map so we can see every city map as we want to.

Second, Call the Google map image

For example, I want to see London Google map. In this case, I can simply use qmap command in ggmap and set the location equal to London.

Then, you can get nice image of London Google map.

But the data I’m using is about crimes in Houston so let’s change it to Houston instead.

Using ‘names’ command, we can get an overview of the variables in the data

For spatial data, as I mentioned in the first paragraph, “lon” and “lat” variables are necessary.

Using ‘dim’ command, we can get the number of rows and columns. Multiplication between rows and columns make dimensions. From Jan 2010 to Aug 2010 in Houston, there were 86,314 crimes. Quite extraordinary!

#Point Data Visualization

Simply, we can use geom_point in ggplot2 package to demonstrate the point map visualization. In this case, I wanted to see the frequencies of different types of crimes.

Pink color is pretty dominant and it indicates that theft is the most predominant crimes in Houston from Jan 2010 to Aug 2010. The second most frequent crime is burglary(the color is confusing, I just hope it’s not murder). Auto theft occurred occasionally.

#Heat Map

If you want to see the density and frequency of the crimes, heat map is the effective.

In this case, we can use stat_density2d for this kind of visualization.

From this heat map, we can observe which area is the most crime-ridden area. Luckily, the campus areas are relatively safer. And it corresponds to the point map that there are lots of points in the first map where it is red in this map. And the red area is the heart of Houston downtown. I hope it has been getting better since then but looks like we’d better be careful around the downtown Houston.

## [R] Harry Potter Sentiment Analysis

Last time, I created word clouds based on Harry Potter. In this post, I will discuss how emotions change throughout each chapter for each book.

This time, you need to download “sentimentr” this time. Lots of useful work can be done by tokenizing at the word level, but sometimes it is useful or necessary to look at different units of text. For example, some sentiment analysis algorithms look beyond only unigrams (i.e. single words) to try to understand the sentiment of a sentence as a whole. These algorithms try to understand that

I’m not having a good day.

is a sad sentence, not a happy one, because of negation. The  `sentimentr` R package are examples of such sentiment analysis algorithms. For these, we may want to tokenize text into sentences.

2. Tokenize text into sentences.

The argument token= sentences attempts to break up text by punctuation.

3. Break up the  text by chapter and sentence.

This will allow us to assess the net sentiment by chapter and by sentence. First, we need to track the sentence numbers and then I create an index that tracks the progress through each chapter. I then unnest the sentences by words. This gives us a tibble that has individual words by sentence within each chapter.

4. Join “afinn” lexicon and compute the net sentiment score

Now, as before, I join the AFINN lexicon and compute the net sentiment score for each chapter.The AFINN lexicon assigns words with a score that runs between -5 and 5, with negative scores indicating negative sentiment and positive scores indicating positive sentiment.

5. Visualize using ggplot

### Result

<Philosopher’s Stone>

This book has the least number of chapters among all seven books. The range of sentiment is from -20 to 15 and it’s the narrowest range of sentiment as well. We can observe that the first chapter is emotionally neutral while chapter 17 contains most emotionally negative and most emotionally positive at the same time. We can see that the ending is relatively happy ending in this book.

<Chamber of Secrets>

It also has the narrowest range of emotions with more chapters. About 25% progress in chapter 1, there is a quite conspicuous negative part and I wonder what it was about.

<Prisoner of Azkaban>

It looks like Prisoner Azkaban does not have many emotionally positive parts. We can see that the highest score is relatively lower than the other two previous series. Instead, the minimum value got lower which indicates that the net sentiment score is lower. Especially, at 50% of chapter 16, we can see dart red color. It indicates that Prisoner of Azkaban got darker than previous ones. But still, it is happy ending.

<Goblet of Fire>

From this book, J.K Rowling started to include more chapters and Goblet of fire has 37 chapters. Emotional range is similar with previous three books. Compared to Prisoner of Azkaban, there are some noticeable blue parts and it may be because Harry Potter getting high scores in Tri Wizard competition was quite exciting. But there are also some red parts which may include Harry Potter being scorned by friends and the death of Cedric. That’s why the ending part is relatively neutral.

<Order of Phoenix>

I feel this one is slightly more colorful than the previous ones.  There are a lot of blues around the middle of the stories but as it goes by, red is pretty dominant. Considering  Sirius Black was killed at the end, it explains why the ending part is not happy ending.

<Half Blood Prince>

It is somehow less colorful than Order of Phoenix. We also should notice that the highest value is the largest in this book. For example, past 75% in chapter 4, the net score is around 30 (I forgot why). Also, there is the darkest red part in chapter 28: 50% and it may be the moment when Dumbledore was killed.

<Deathly Hallows>

Interestingly, this book has the lowest net score: -40. In Chapter 17 after 50% progression, there is -40 part. According to the story, it is the part where Harry confronted Bathilda changing into a snake. Besides that, we can see that the negative and neutral sentiment is dominant in this one. But we know that it ends well!

## [Python]Principal Component Analysis and K-means clustering with IMDB movie datasets

Hello, today’s post would be the first post that I present the result in Python! Although I love R and I’m loyal to it, Python is widely loved by many data scientists.  Python is quite easy to learn and it has a lot of great functions.

In this post, I implemented unsupervised learning methods: 1. Principal Component Analysis and 2. K-means Clustering. Then a reader who has no background knowledge in Machine Learning would think,”what the hell is unsupervised learning?” I will try my best to explain this concept

### Unsupervised Learning

Ok, let’s imagine you are going to backpacking to a new country. Isn’t it exciting? But you did not know much about the country – their food, culture, language etc. However from day 1, you start making sense there, learning to eat new cuisines including what not to eat, find a way to that beach.

In this example,you have lots of information but you do not know what to do with it initially. There is no clear guidance and you have to find the way by yourself. Like this traveling example, unsupervised learning is the method of training your machine learning task only with a set of inputs. Principal Component Analysis and K-means clustering are the most famous examples of unsupervised learning. I will explain them a little bit later.

### Data

Before I begin talking about how I analyzed the data, let’s talk about the data. There are total 5,043 movies with 28 attributes. The attributes range from director name to the number of facebook likes.

### 1. Data Cleaning

In Statistics class, we often get clean data: no missing values, no NA values. But in reality, the clean data is just like a dream. There are always some messed part of the data and it’s our job to trim the data useable before executing the analysis.

Here are some libraries you need for this post.

First, let’s do some filtering to extract only the numbered columns and not the ones with words. So, I created a Python list containing the numbered column names “num_list”

By the way, when it comes to using Python, pandas library is a must-have item. Using pandas library, we can create a new dataframe (movie_num) containing just the numbers

By using function “fillna(filtering NA)”, we can easily discard NaN values.

If the distribution of certain variables are skewed, we can implement standardization.

### 2. Correlation Analysis

#### Hexbin Plot

Let’s look at some hexbin visualisations first to get a feel for how the correlations between the different features compare to one another. In the hexbin plots, the lighter in color the hexagonal pixels, the more correlated one feature is to another.

This is a Hexbin Plot between IMDB Scroe and gross revenue. We can see it’s lighter around the score between 6 and 7.

This is a Hexbin Plot between IMDB Scroe and duration(days). Again, the score between 6 and 7 is lighter.

We can examine the correlation more using Pearson correlation plot.

As we can see from the heatmap, there are regions (features) where we can see quite positive linear correlations amongst each other, given the darker shade of the colours – top left-hand corner and bottom right quarter. This is a good sign as it means we may be able to find linearly correlated features for which we can perform PCA projections on.

### 3. EXPLAINED VARIANCE MEASURE &Principal Component Analysis

Now you know what unsupervised learning is (I hope so). Then, let me explain about principal component analysis. The explanation would not be as entertaining as the one in unsupervised learning but I’ll try my best!

Principal component analysis (PCA) is a technique used to emphasize variation and bring out strong patterns in a dataset. It’s often used to make data easy to explore and visualize.  Principal components are dimensions along which your data points are most spread out:

<From: https://www.quora.com/What-is-an-intuitive-explanation-for-PCA>

Let me give you an example. Imagine that you are a nutritionist trying to explore the nutritional content of food. What is the best way to differentiate food items? By vitamin content? Protein levels? Or perhaps a combination of both?

Knowing the variables that best differentiate your items has several uses:

1. Visualization. Using the right variables to plot items will give more insights.

2. Uncovering Clusters. With good visualizations, hidden categories or clusters could be identified. Among food items for instance, we may identify broad categories like meat and vegetables, as well as sub-categories such as types of vegetables.

The question is, how do we derive the variables that best differentiate items?

So, the first step to answer this question is Principal Component Analysis.

A principal component can be expressed by one or more existing variables. For example, we may use a single variable – vitamin C – to differentiate food items. Because vitamin C is present in vegetables but absent in meat, the resulting plot (below, left) will differentiate vegetables from meat, but meat items will clumped together.

To spread the meat items out, we can use fat content in addition to vitamin C levels, since fat is present in meat but absent in vegetables. However, fat and vitamin C levels are measured in different units. So to combine the two variables, we first have to normalize them, meaning to shift them onto a uniform standard scale, which would allow us to calculate a new variable – vitamin C – fat. Combining the two variables helps to spread out both vegetable and meat items.

The spread can be further improved by adding fiber, of which vegetable items have varying levels. This new variable – (vitamin C + fiber) – fat – achieves the best data spread yet.

So,  that’s my explanation of Principal Component analysis and K-means clustering at the same time. Let me apply Principal Component Analysis to this dataset and show how it works.

#### Explained Variance Measure

I will be using a particular measure called Explained Variance which will be useful in this context to help us determine the number of PCA projection components we should be looking at.

Before calculating explained variance, we need to get eigenvectors and eigenvalues.The eigenvectors and eigenvalues of a covariance (or correlation) matrix represent the “core” of a PCA: The eigenvectors (principal components) determine the directions of the new feature space, and the eigenvalues determine their magnitude. In other words, the eigenvalues explain the variance of the data along the new feature axes.

After sorting the eigenpairs, the next question is “how many principal components are we going to choose for our new feature subspace?”. The explained variance tells us how much information (variance) can be attributed to each of the principal components.

From the plot above, it can be seen that approximately 90% of the variance can be explained with the 9 principal components. Therefore for the purposes of this notebook, let’s implement PCA with 9 components ( although to ensure that we are not excluding useful information, one should really go for 95% or greater variance level which corresponds to about 12 components).

There does not seem to be any discernible clusters. However keeping in mind that our PCA projections contain another 7 components, perhaps looking at plots with the other components may be fruitful. For now, let us assume that will be trying a 3-cluster (just as a naive guess) KMeans to see if we are able to visualize any distinct clusters.

### 5.Visualization with K-means clustering

This KMeans plot looks more promising now as if our simple clustering model assumption turns out to be right, we can observe 3 distinguishable clusters via this color visualization scheme. However I would also like to generate a KMeans visualization for other possible combinations of the projections against one another. I will use Seaborn’s convenient pairplot function to do the job. Basically pairplot automatically plots all the features in the dataframe (in this case our PCA projected movie data) in pairwise manner. I will pairplot the first 3 projections against one another and the resultant plot is given below: