Thursday, 5 March 2015

Part 3a: Plotting with ggplot2


We will start off this first section of Part 3 with a brief introduction of the plotting system ggplot2. Then, with the attention focused mainly on the syntax, we will create a few graphs, based on the weather data we have prepared previously.

Next, in Part 3b, where we will be doing actual EDA, specific visualisations using ggplot2 will be developed, in order to address the following question: Are there any good predictors, in our data, for the occurrence of rain on a given day?

Lastly, in Part 4, we will use Machine Learning algorithms to check to which extent rain can be predicted based on the variables available to us.

But for now, let's see what this ggplot2 system is all about and create a few nice looking figures.

 

The concept of ggplot2


The R package ggplot2, created by Hadley Wickham, is an implementation of Leland Wilkinson's Grammar of Graphics, which is a systematic approach to describe the components of a graphic. In maintenance mode (i.e., no active development) since February 2014, ggplot2 it is the most downloaded R package of all time.

In ggplot2, the graphics are constructed out of layers. The data being plotted, coordinate system, scales, facets, labels and annotations, are all examples of layers. This means that any graphic, regardless of how complex it might be, can be built from scratch by adding simple layers, one after another, until one is satisfied with the result. Each layer of the plot may have several different components, such as: data, aesthetics (mapping), statistics (stats), etc. Besides the data (i.e., the R data frame from where we want to plot something), aesthetics and geometries are particularly important, and we should make a clear distinction between them:
  • Aesthetics are visual elements mapped to data; the most common aesthetic attributes are the x and y values, colour, size, and shape;
  • Geometries are the actual elements used to plot the data, for example: point, line, bar, map, etc.
So, what does this mean in practice? Suppose we want to draw a scatter plot to show the relationship between the minimum and maximum temperatures on each day, controlling for the season of the year using colour (good for categorical variables), and the rain using size (good for continuous variables). We would map the aesthetics x, y, colour and size to the respective data variables, and then draw the actual graphic using a geometry - in this case, geom_point().

Here are some important aspects to have in mind:
  • Not all aesthetics of a certain geometry have to be mapped by the user, only the required ones. For example, in the case above we didn't map the shape of the point, and therefore the system would use the default value (circle); it should be obvious, however, that x and y have no defaults and therefore need to be either mapped or set;
  • Some aesthetics are specific of certain geometries - for instance, ymin and ymax are required aesthetics for geom_errorbar(), but don't make sense in the context of geom_point();
  • We can map an aesthetic to data (for example, colour = season), but we can also set it to a constant value (for example, colour = "blue"). There are a few subtleties when writing the actual code to map or set an aesthetic, as we will see below.
 

 

The syntax of ggplot2


A basic graphic in ggplot2 consists of initializing an object with the function ggplot() and then add the geometry layer. 
  • ggplot() - this function is always used first, and initializes a ggplot object. Here we should declare all the components that are intended to be common to all the subsequent layers. In general, these components are the data (the R data frame) and some of the aesthetics (mapping visual elements to data, via the aes() function). As an example:
ggplot(data = dataframe, aes(x = var1, y= var2, colour = var3, size = var4))
  • geom_xxx() - this layer is added after ggplot() to actually draw the graphic If all the required aesthetics were already declared in ggplot(), it can be called without any arguments. If an aesthetic previously declared in ggplot() is also declared in the geom_xxx layer, the latter overwrites the former. This is also the place to map or set specific components specific to the layer. A few real examples from our weather data should make it clear:
# This draws a scatter plot where season controls the colour of the point
ggplot(data = weather, aes(x = l.temp, y= h.temp, colour = season) + geom_point()

# This draws a scatter plot but colour is now controlled by dir.wind (overwrites initial definition)
ggplot(data = weather, aes(x = l.temp, y= h.temp, colour = season) + geom_point(aes(colour=dir.wind))

# This sets the parameter colour to a constant instead of mapping it to data. This is outside aes()!
ggplot(data = weather, aes(x = h.temp, y= l.temp) + geom_point(colour = "blue")

# Mapping to data (aes) and setting to a constant inside the geom layer
ggplot(data = weather, aes(x = l.temp, y= h.temp) + geom_point(aes(size=rain), colour = "blue")

# This is a MISTAKE - mapping with aes() when intending to set an aesthetic. It will create a column in the data named "blue", which is not what we want.
ggplot(data = weather, aes(x = l.temp, y= h.temp) + geom_point(aes (colour = "blue"))

When learning the ggplot2 system, this is a probably a good way to make the process easier:
  1. Decide what geom you need - go here, identify the geom best suited for your data, and then check which aesthetics it understands (both the required, that you will need to specify, and the optional ones, that more often than not are useful);
  2. Initialize the plot with ggplot() and map the aesthetics that you intend to be common to all the subsequent layers (don't worry too much about this, you can always overwrite aesthetics if needed);
  3. Make sure you understand if you want an aesthetic to be mapped to data or set to a constant. Only in the former the aes() function is used. 
Although aesthetics and geometries are probably the most important part of ggplot2, there are more layers/components that we will be using often, such as scales (both colour scales and axes-controlling ones), statistics, and faceting (this one surely in Part 3b). Now, let's make some plots using different geometries!


 

Plotting a time series - point geom with smoother curve

 
 
# Time series of average daily temperature, with smoother curve

ggplot(weather,aes(x = date,y = ave.temp)) +
  geom_point(colour = "blue") +
  geom_smooth(colour = "red",size = 1) +
  scale_y_continuous(limits = c(5,30), breaks = seq(5,30,5)) +
  ggtitle ("Daily average temperature") +
  xlab("Date") +  ylab ("Average Temperature ( ºC )")
 

 

Instead of setting the colour to blue, we can map it to the value of the temperature itself, using the aes() function. Moreover, we can define the color gradient we want - in this case, I told that cold days should be blue and warm days red. I forced the gradient to go through green, and pure green should occur when the temperature is 16 ºC. There are many ways to work with colour, and it would be impossible for us to cover everything. 

 
# Same but with colour varying

ggplot(weather,aes(x = date,y = ave.temp)) + 
  geom_point(aes(colour = ave.temp)) +
  scale_colour_gradient2(low = "blue", mid = "green" , high = "red", midpoint = 16) + 
  geom_smooth(color = "red",size = 1) +
  scale_y_continuous(limits = c(5,30), breaks = seq(5,30,5)) +
  ggtitle ("Daily average temperature") +
  xlab("Date") +  ylab ("Average Temperature ( ºC )")
 
 
The smoother curve (technically a loess) drawn on the graphs above shows the typical pattern for a city in the northern hemisphere, i.e., higher temperatures in July (summer) and lower in January (winter).

 

 

Analysing the temperature by season - density geom

 
 
# Distribution of the average temperature by season - density plot

ggplot(weather,aes(x = ave.temp, colour = season)) +
  geom_density() +
  scale_x_continuous(limits = c(5,30), breaks = seq(5,30,5)) +
  ggtitle ("Temperature distribution by season") +
  xlab("Average temperature ( ºC )") +  ylab ("Probability")

 


I find this graph interesting. Spring and autumn seasons are often seen as the transition from cold to warm and warm to cold days, respectively. The spread of their distributions reflect the high thermal amplitude of these seasons. On the other hand, winter and summer average temperatures are much more concentrated around a few values, and hence the peaks shown on the graph.

 

Analysing the temperature by month - violin geom with jittered points overlaid


# Label the months - Jan...Dec is better than 1...12

weather$month = factor(weather$month,
                       labels = c("Jan","Fev","Mar","Apr",
                                  "May","Jun","Jul","Aug","Sep",
                                  "Oct","Nov","Dec"))

# Distribution of the average temperature by month - violin plot,
# with a jittered point layer on top, and with size mapped to amount of rain

ggplot(weather,aes(x = month, y = ave.temp)) +
  geom_violin(fill = "orange") +
  geom_point(aes(size = rain), colour = "blue", position = "jitter") +
  ggtitle ("Temperature distribution by month") +
  xlab("Month") +  ylab ("Average temperature ( ºC )")
 
 
A violin plot is sort of a mixture of a box plot with a histogram (or even better, a rotated density plot). I also added a point layer on top (with jitter for better visualisation), in which the size of the circle is mapped to the continuous variable representing the amount of rain. As can be seen, there were several days in the autumn and winter of 2014 with over 60 mm of precipitation, which is more than enough to cause floods in some vulnerable zones. In subsequent parts, we will try to learn more about rain and its predictors, both with visualisation techniques and machine learning algorithms.
 

 

Analysing  the correlation between low and high temperatures

 
 
# Scatter plot of low vs high daily temperatures, with a smoother curve for each season

ggplot(weather,aes(x = l.temp, y = h.temp)) +
  geom_point(colour = "firebrick", alpha = 0.3) + 
  geom_smooth(aes(colour = season),se= F, size = 1.1) +
  ggtitle ("Daily low and high temperatures") +
  xlab("Daily low temperature ( ºC )") +  ylab ("Daily high temperature ( ºC )") 



The scatter plot shows a positive correlation between the low and high temperatures for a given day, and it holds true regardless of the season of the year. This is a good example of a graphic where we set an aesthetic for one geom (point colour) and mapped an aesthetic for another (smoother curve colour).

 

Distribution of low and high temperatures by the time of day

 

In this final example, we would like to plot the frequencies of two data series (i.e. two variables) on the same graph, namely, the lowest and highest daily temperatures against the hour when they occurred. The way these variables are represented in our data set is called wide data form, which is used, for instance, in MS Excel, but it is not the first choice in ggplot2, where the long data form is preferred.

Let's think a bit about this important, albeit theoretical, concept. When we plotted the violin graph before, did we have wide data (12 numerical variables, each one containing the temperatures for each month), or long data (1 numerical variable with all temperatures and 1 grouping variable with the month information)? Definitely the latter. This is exactly need what we need to do to our 2 temperature variables: create a grouping variable that, for each row, represents either a low temperature or a high temperature, and a numerical variable indicating the corresponding hour of the day. The example in the code below should help understanding this concept (for further information on this topic, please have a look at the tidy data article mentioned in Part 1 of this tutorial).

> # We will use the melt function from reshape2 to convert from wide to long form
> library(reshape2) 
 
> # select only the variables that are needed (day and temperatures) and assign to a new data frame
> temperatures <- weather[c("day.count","h.temp.hour","l.temp.hour")] 
 
> # Temperature variables in columns
> head(temperatures)
  day.count h.temp.hour l.temp.hour
1         1           0           1
2         2          11           8
3         3          14          21
4         4           2          11
5         5          13           2
6         6           0          20 
 
> dim(temperatures)
[1] 365   3 
 
> # The temperatures are molten into a single variable called l.h.temp 
> temperatures <- melt(temperatures,id.vars = "day.count",
                       variable.name = "l.h.temp", value.name = "hour")

> # See the difference?
> head(temperatures)
  day.count    l.h.temp hour
1         1 h.temp.hour    0
2         2 h.temp.hour   11
3         3 h.temp.hour   14
4         4 h.temp.hour    2
5         5 h.temp.hour   13
6         6 h.temp.hour    0 
 
> tail(temperatures)
    day.count    l.h.temp hour
725       360 l.temp.hour    7
726       361 l.temp.hour    4
727       362 l.temp.hour   23
728       363 l.temp.hour    8
729       364 l.temp.hour    5
730       365 l.temp.hour    8 
 
> # We now have twice the rows. Each day has an observation for low and high temperatures
> dim(temperatures)
[1] 730   3 
 
> # Needed to force the order of the factor's level
> temperatures$hour <- factor(temperatures$hour,levels=0:23)

# Now we can just fill the colour by the grouping variable to visualise the two distributions 
ggplot(temperatures) +
  geom_bar(aes(x = hour, fill = l.h.temp)) +
  scale_fill_discrete(name= "", labels = c("Daily high","Daily low")) +
  scale_y_continuous(limits = c(0,100)) +
  ggtitle ("Low and high temperatures - time of the day") +
  xlab("Hour") +  ylab ("Frequency")
 
 

The graph clearly shows with we know from experience: the daily lowest temperatures tend to occur early in the morning and the highest in the afternoon. But can you see those outliers in the distribution? There were a few days were the highest temperature was during the night, and a few others where the lowest was during the day. Can this somehow be correlated to rain? This is the kind of questions we will be exploring in Part 3b, where we will use EDA techniques (visualisations) to identify potential predictors for the occurrence of rain. Then, in Part 4, using data mining and machine learning algorithms, not only we will find the best predictors for rain (if any) and determine the accuracy of the models, but also the individual contribution of each variable in our data set.

More to come soon... stay tuned!

Friday, 20 February 2015

Part 2: Data Preparation


In Part 1 I have introduced the weather data set we will be using in this series of tutorials. We are now going to have the data prepared for the subsequent EDA. We will recode and transform variables, change their types, and perform some basic data checks. Feel free to follow along with the analysis (click here to download the weather data), bearing in mind you can type ?function_name to get help about some specific R function, for instance, ?head

Importing the data


To start off, let's read in the data to an R data frame and run some basic commands:

# Make sure the file is in your current working directory
 
> weather <- read.csv("weather_2014.csv",sep=";",stringsAsFactors=FALSE)
 
> dim(weather)
[1] 365  14

> names(weather)
 [1] "day.count"      "day"            "month"          "season"        
 [5] "l.temp"         "h.temp"         "ave.temp"       "l.temp.time"   
 [9] "h.temp.time"    "rain"           "ave.wind"       "gust.wind"     
[13] "gust.wind.time" "dir.wind"
 
> head(weather)
  day.count day month season l.temp h.temp ave.temp l.temp.time h.temp.time rain
1         1   1     1 Winter   12.7   14.0     13.4       01:25       23:50 32.0
2         2   2     1 Winter   11.3   14.7     13.5       07:30       11:15 64.8
3         3   3     1 Winter   12.6   14.7     13.6       21:00       14:00 12.7
4         4   4     1 Winter    7.7   13.9     11.3       10:35       01:50 20.1
5         5   5     1 Winter    8.8   14.6     13.0       01:40       12:55  9.4
6         6   6     1 Winter   11.8   14.4     13.1       19:35       00:05 38.9
  ave.wind gust.wind gust.wind.time dir.wind
1     11.4      53.1          15:45        S
2      5.6      41.8          22:25        S
3      4.3      38.6          00:00      SSW
4     10.3      66.0          09:05       SW
5     11.6      51.5          13:50      SSE
6      9.9      57.9          08:10      SSE 
 
 
It seems we have correctly loaded the data into R. Notice that the variables in the data file are separated not by a comma, but by a semicolon, and hence the need to set the sep = ";" argument. We also told R not to import strings as factors (i.e., categorical variables). In many cases some character variables are indeed strings and some integer variables are actually factors. This means there is almost always manual work to be done after importing, and therefore I prefer to read in the variables without any initial processing.

Let's now have a look at the structure of the weather data frame, using one of the most useful R functions, str().

> str(weather)
'data.frame': 365 obs. of  14 variables:
 $ day.count     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ day           : int  1 2 3 4 5 6 7 8 9 10 ...
 $ month         : int  1 1 1 1 1 1 1 1 1 1 ...
 $ season        : chr  "Winter" "Winter" "Winter" "Winter" ...
 $ l.temp        : num  12.7 11.3 12.6 7.7 8.8 11.8 11.4 12.4 9.2 8.3 ...
 $ h.temp        : num  14 14.7 14.7 13.9 14.6 14.4 14.8 15.6 18.4 14.8 ...
 $ ave.temp      : num  13.4 13.5 13.6 11.3 13 13.1 13.5 14.1 12.9 11 ...
 $ l.temp.time   : chr  "01:25" "07:30" "21:00" "10:35" ...
 $ h.temp.time   : chr  "23:50" "11:15" "14:00" "01:50" ...
 $ rain          : num  32 64.8 12.7 20.1 9.4 38.9 2 1.5 0 0 ...
 $ ave.wind      : num  11.4 5.6 4.3 10.3 11.6 9.9 6.6 5.9 0.2 1.4 ...
 $ gust.wind     : num  53.1 41.8 38.6 66 51.5 57.9 38.6 33.8 16.1 24.1 ...
 $ gust.wind.time: chr  "15:45" "22:25" "00:00" "09:05" ...
 $ dir.wind      : chr  "S" "S" "SSW" "SW" ...

Based on the output of this function, and having in the mind the goal of producing visualisations and potentially build models, how would you change the variables in the data set? Here are my thoughts:
  • Day and month are coded as integers but they should be factors (categorical variables); the same applies to the character variables representing the season and wind direction;
  • Looking at the first values for the wind direction - "S" "S" "SSW" "SW" - it seems that a 16-wind compass rose has been used. Since there are only 365 days in the year, do we have sufficient observations for each of the 16 groups, or could we try to group them into 8 principal winds or even only the 4 cardinal directions?
  • The day count (number of days since the beginning of the year) is useful, but we would like to have dates on the x axis when plotting instead of an index. This variable should therefore be transformed;
  • The three time variables show the exact minute where the corresponding event occurred. We would most likely benefit by doing some aggregation, and rounding to the nearest hour seems a good option. After that we would convert the hour variable to factor.

It is also quite common to check for missing values (coded as NA by default). As with almost everything, there are many ways to do it in R. Here are just two of them:

# One way (sum of NA over the entire data set)
> sum(is.na(weather))
[1] 0

# Another way (number of complete observations)
> nrow(weather)
[1] 365 
 
> sum(complete.cases(weather))
[1] 365 
 
> nrow(weather) == sum(complete.cases(weather))
[1] TRUE
 

 

Create factors

 

Variables can be easily coerced to factors using the as.factor() function, which is an abbreviated form of the main function, factor(). The first one, however, will order the levels by the alphabet when the original variable is a string (class character in R), which might be something we don't want for some variables, for example, the season of the year. The code below shows how factors are created and why they differ from the other types of variables.

> # Before (365 independent strings)
> class(weather$season)
[1] "character"

> summary(weather$season)
   Length     Class      Mode 
      365 character character 

> weather$season <- factor(weather$season,
                    levels = c("Spring","Summer","Autumn","Winter"))

> # After (4 categories, ordered by "levels")
> class(weather$season)
[1] "factor"

> summary(weather$season)
Spring Summer Autumn Winter 
    92     92     91     90
 
> # Using as.factor() when the order doesn't matter or original var. is integer

> weather$day <- as.factor(weather$day)
> weather$month <- as.factor(weather$month)
> weather$dir.wind <- as.factor(weather$dir.wind)
 

 

Dealing with the wind


Let's start by checking whether there are actually 16 directions in the dir.wind variable and, if so, determine whether we have sufficient number of observations in each group.

> # Number of unique values
> length(unique(weather$dir.wind))
[1] 16 
 
> # Absolute frequency (table function)
> table(weather$dir.wind)

  E ENE ESE   N  NE NNE NNW  NW   S  SE SSE SSW  SW   W WNW WSW 
 11  15   2  18  25   8  37 108  26  24  31  17  11   5  24   3 
 
> # Making it relative (prop.table function)
> rel <- round(prop.table(table(weather$dir.wind))*100,1)
> rel

   E  ENE  ESE    N   NE  NNE  NNW   NW    S   SE  SSE  SSW   SW    W  WNW  WSW 
 3.0  4.1  0.5  4.9  6.8  2.2 10.1 29.6  7.1  6.6  8.5  4.7  3.0  1.4  6.6  0.8 
 
> # Bringing some order to the table
> sort(rel,decreasing = TRUE)

  NW  NNW  SSE    S   NE   SE  WNW    N  SSW  ENE    E   SW  NNE    W  WSW  ESE 
29.6 10.1  8.5  7.1  6.8  6.6  6.6  4.9  4.7  4.1  3.0  3.0  2.2  1.4  0.8  0.5
 
It can be seen that the relative frequency is less than 5% for more than half of the groups. Unfortunately, we don't have the actual wind direction in degrees. But, making use of some domain knowledge and the information in the table above, let's try to give a reasonable answer to the following question: Is it more likely for a value in the NNW group to be closer to NW to N?
(assuming the direction isn't exactly the midpoint between NW and N, in which case it would be a pure NNW). The numbers show that NW would be far more likely. Even though this logic doesn't necessarily apply to all of the directions we are trying to eliminate, let's assume this criterion of recode them as the closest ordinal direction (by definition, "NW","NE","SE", "SW" are called ordinal) in a new variable.

As long as the analyst knows to explain why and how some new variable was created, and bears in mind it may lack accuracy, it is perfectly fine to add it to the data set. It may be useful or useless, and that is something he will try to figure out during the stages of visualisation or modelling.

> # Transforming wind direction variable: from 16 to 8 principal winds 
 
> # Create a copy from the original variable...
> weather$dir.wind.8 <- weather$dir.wind 
 
> # ...and then simply recode some of the variables
> weather$dir.wind.8 <- ifelse(weather$dir.wind %in%  c("NNE","ENE"),
                               "NE",as.character(weather$dir.wind.8)) 
 
> weather$dir.wind.8 <- ifelse(weather$dir.wind %in% c("NNW","WNW"),
                               "NW",as.character(weather$dir.wind.8)) 
 
> weather$dir.wind.8 <- ifelse(weather$dir.wind %in% c("WSW","SSW"),
                               "SW",as.character(weather$dir.wind.8)) 
 
> weather$dir.wind.8 <- ifelse(weather$dir.wind %in% c("ESE","SSE"),
                               "SE",as.character(weather$dir.wind.8)) 
 
> # create factors, ordered by "levels" 
> weather$dir.wind.8 <- factor(weather$dir.wind.8,
                        levels = c("N","NE","E","SE","S","SW","W","NW"))
 
 
> # Checking the length of the new variable 
> length(unique(weather$dir.wind.8))
[1] 8 
 
> # A 2-way table (direction vs season), with relative frequencies calculated
    over margin = 2 (i.e., the columns)  
 
> round(prop.table(table(weather$dir.wind.8,weather$season),margin = 2)*100,1)
    
     Spring Summer Autumn Winter
  N     1.1    3.3   12.1    3.3
  NE   14.1    5.4   20.9   12.2
  E     0.0    0.0    5.5    6.7
  SE   13.0   14.1   20.9   14.4
  S     5.4   12.0    4.4    6.7
  SW    6.5    8.7    2.2   16.7
  W     2.2    0.0    1.1    2.2
  NW   57.6   56.5   33.0   37.8 


The function ifelse() is one of the classical ways to populate a new variable based on the value, or calculation over the value, of any other(s). When the original variable is numeric, cut() is often simpler and used instead.

Just as a side note, it is uncommon and considered bad practice to use for loops and if statements in R when the goal is to loop through the rows of a column and apply some function when a certain condition is met. R supports vectorisation, which is a much more efficient way to accomplish the same thing. In fact, the ifelse() function is the vectorised way that makes the for-if-else construct unnecessary.

 

 

We need a date


To create a date in R, we just need to pass to the function as.Date() a string with an appropriate format. We can then add and subtract days using the usual math operators. Here is the code to calculate the date based on the day.count variable.

> first.day <- "2014-01-01"
> class(first.day)
[1] "character" 
 
> first.day <- as.Date(first.day)
> class(first.day)
[1] "Date" 
 
> Here is where we actually calculate the date 
> weather$date  <- first.day + weather$day.count - 1 
 
> head(weather$day.count)
[1] 1 2 3 4 5 6 
 
> head(weather$date)
[1] "2014-01-01" "2014-01-02" "2014-01-03" "2014-01-04" "2014-01-05" "2014-01-06"

 

 

Not so hard times


The last thing we need to do is to round (to the nearest hour) the time at which a certain event occurred  (lower temperature, higher temperature, and wind gust). Working with times in R is a bit more complicated than working with dates, with the former having two alternative classes to represent it: POSIXct and POSIXlt. The first stores the date and time as a simple number, representing the seconds since the UNIX epoch (Jan 1, 1970); the second stores the date and time in a list, with elements for seconds, hours, years, among others. Since we are interested in extracting the hour information, after rounding, we will use the more complete POSIXlt class.

> # Store date and time as POSIXlt class
> l.temp.time.date <- as.POSIXlt(paste(weather$date,weather$l.temp.time))
> head(l.temp.time.date)
[1] "2014-01-01 01:25:00 GMT" "2014-01-02 07:30:00 GMT" "2014-01-03 21:00:00 GMT"
[4] "2014-01-04 10:35:00 GMT" "2014-01-05 01:40:00 GMT" "2014-01-06 19:35:00 GMT" 
 
> # Round to the nearest hour
> l.temp.time.date <- round(l.temp.time.date,"hours")
> head(l.temp.time.date)
[1] "2014-01-01 01:00:00 GMT" "2014-01-02 08:00:00 GMT" "2014-01-03 21:00:00 GMT"
[4] "2014-01-04 11:00:00 GMT" "2014-01-05 02:00:00 GMT" "2014-01-06 20:00:00 GMT" 
 
> # Which attributes are stored in the POSIXlt time variable?
> attributes(l.temp.time.date)
$names
 [1] "sec"    "min"    "hour"   "mday"   "mon"    "year"   "wday"   "yday"   "isdst" 
[10] "zone"   "gmtoff"

$class
[1] "POSIXlt" "POSIXt" 

$tzone
[1] ""    "GMT" "BST"

> # Extract the value of the hour attribute as a number and add it to the data set
> weather$l.temp.hour <- l.temp.time.date [["hour"]]
> head(weather$l.temp.hour)
[1]  1  8 21 11  2 20 
 
> # Lastly, the integer is converted to factor
> weather$l.temp.hour <- as.factor(weather$l.temp.hour)
> head(weather$l.temp.hour)
[1] 1  8  21 11 2  20
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 17 18 19 20 21 22 23

 

 

The prepared data set


After all the processing we have done, let's call str() again to see what our final data set looks like. We now have a date variable (to plot a time series) and several factor variables (commonly used to identify different groups of a numerical variable when creating a data visualisation).

> str(weather)
'data.frame': 365 obs. of  19 variables:
 $ day.count     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ day           : Factor w/ 31 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ month         : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ season        : Factor w/ 4 levels "Spring","Summer",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ l.temp        : num  12.7 11.3 12.6 7.7 8.8 11.8 11.4 12.4 9.2 8.3 ...
 $ h.temp        : num  14 14.7 14.7 13.9 14.6 14.4 14.8 15.6 18.4 14.8 ...
 $ ave.temp      : num  13.4 13.5 13.6 11.3 13 13.1 13.5 14.1 12.9 11 ...
 $ l.temp.time   : chr  "01:25" "07:30" "21:00" "10:35" ...
 $ h.temp.time   : chr  "23:50" "11:15" "14:00" "01:50" ...
 $ rain          : num  32 64.8 12.7 20.1 9.4 38.9 2 1.5 0 0 ...
 $ ave.wind      : num  11.4 5.6 4.3 10.3 11.6 9.9 6.6 5.9 0.2 1.4 ...
 $ gust.wind     : num  53.1 41.8 38.6 66 51.5 57.9 38.6 33.8 16.1 24.1 ...
 $ gust.wind.time: chr  "15:45" "22:25" "00:00" "09:05" ...
 $ dir.wind      : Factor w/ 16 levels "E","ENE","ESE",..: 9 9 12 13 11 11 10 10 4 7 ...
 $ dir.wind.8    : Factor w/ 8 levels "N","NE","E","SE",..: 5 5 6 6 4 4 4 4 1 8 ...
 $ date          : Date, format: "2014-01-01" "2014-01-02" ...
 $ l.temp.hour   : Factor w/ 20 levels "0","1","2","3",..: 2 9 18 12 3 17 8 1 8 9 ...
 $ h.temp.hour   : Factor w/ 19 levels "0","1","2","3",..: 1 8 11 3 10 1 12 10 11 9 ...
 $ gust.wind.hour: Factor w/ 24 levels "0","1","2","3",..: 17 23 1 10 15 9 13 15 15 15 ...

  

 

Final notes on data "wrangling" and R


As we have seen in Part 1, the process of cleaning and transforming the raw data is almost always required prior to starting the actual analysis. In fact, in many real life cases, the visualisation and modeling stages are easier and less time-consuming than having the data ready to be explored.

The R language is extremely powerful to create visualisations and build models. It is often considered, however, that the language is not the most user-friendly when it comes to prepare the data (still true, but not as much as a few years ago). Here are some of the alternatives we, the analysts, have at our disposal:
  • For reasonably simple and tidy data sets, like the one we have been using in this tutorial, Excel is usually sufficient; a combination of Pivot Tables, Vlookup() and/or Index()/Match() and a few basic formatting functions would have accomplished the same we have done here;
  • When some more advanced processing is required (for instance, the use of regex), and even though R supports regular expressions and provides functions in its base package, some prefer to use other languages (in this case, Perl or Python would be good options);
  • An ever increasing alternative option is to use other R packages that provide wrapper functions for the ones in its base. For instance, we would probably have dealt with our three time variables easily using the lubridate package than using the base functions; had we needed to format strings and the stringr package would provide us with several consistent wrappers, making it simpler to process the text.

Onward and upward to Part 3 of this series of tutorials, where we will be creating a few  visualisations to gain insight from our weather data.        

Saturday, 14 February 2015

Part 1: Introduction


The ultimate goal of every data scientist is to extract as much valuable information as possible from a given data set. We want to be able to predict the future based on the past, to discover very deep and hidden patterns in the data, and to expand the current base of knowledge in some specific domain. With this in mind, several machine learning algorithms, from Neural Networks to Support Vector Machines, from Naïve Bayes to Random Forests, have been developed over the years. In many situations, when correctly applied, these can provide greater insight from the data than any human could do, how clever they might be, based on his analytical skills.
 
Daily experience shows, however, that the process of analysing data seldom goes past the Exploratory Data Analysis (EDA). In fact, more often than not, the classical EDA approach of summarising and plotting is sufficient for the analyst to build a solid intuition of what the data is trying to convey, and sometimes raising new questions that might be addressed afterwards, if a model is to be developed.

Unfortunately, in many cases we cannot perform EDA right away. Raw data is often messy, unstructured, badly coded, inconsistent, or just plain wrong. It is often said that the analyst spends 80% of the time preparing the data, and only 20% actually doing analysis and modelling. This initial data wrangling consists of handling missing values, remove duplicates, transforming variables, format variable types, recode values, detect outliers to ascertain data integrity, and others. Our goal is to have what is called tidy data by the end of the process.

Data analysis using R


In this series of intermediate-level tutorials, I will guide you through the process of analysing an actual data set. We will start by preparing the data and then use a few EDA techniques to get a grasp of what’s in it. To accomplish this, we will be using the lingua franca of data science, the R programming language. There are many advantages and disadvantages to using R, and I summarise below the ones I personally deem the most relevant.

Pros
  • There is a big community of R users – as of 2014, there were more than 5 000 user-contributed packages only in the main repository (CRAN), and around 150 000 R functions (software popularity). If you need a function for some specific purpose, it is highly likely that someone else has already created it before. If you have any doubt about the R syntax, you’ll probably find the answer online easily (stackoverflow);
  • The superb graphics capabilities offered by the ggplot2 package - this plotting system implements the grammar of graphics, a new way of thinking about the visual representation of the data. If you have to pick a single topic lo learn in R, ggplot2 is arguably the best option. It is a language by itself that will allow you to produce high quality plots in a short amount of time;
  •  It’s free! (This is not unique to R, though).

Cons
  • R has a steep learning curve – the paradigm is different from the mainstream languages and even from other statistical packages. It is highly interactive, where to complete an analysis you often call one function, take the result and use it to feed the next function and so on, in a cycle that can be quite extensive;
  • Some of the syntax is far from intuitive and even cumbersome – for example, using the intuitive sort() function to sort a data frame yields in nasty results; you need to use the order() function instead, not directly, but in a rather convoluted way. It must be said, however, that several packages have been developed to make the analysis much easier, namely plyr/dplyr to manipulate (filter, transform, summarise) data, lubridate and stringr to lessen the burden when dealing with dates and strings, respectively, sqldf to run SQL statements directly over R data frames, among a few other packages;
  • R can be a bit slow, but that’s more in the context of developing algorithms than about performing interactive EDA. In fact, the bottleneck of the process is the analyst himself. The time we spend thinking about what information to extract and how we want to visualise it is several orders of magnitude greater than the time it takes R to actually draw the graphics.

Without further ado, let’s have a look at the data set we will be using in this series of tutorials.

The weather data set


The data set consists of daily records of several meteorological parameters, measured in the city of Porto over the year of 2014. We have, then, 365 observations for each of the following 14 variables:

day.count – number of days passed since the beginning of the year
day – day of the month
month – month of the year
season – season of the year
l.temp, h.temp, ave.temp – lowest, highest and average temperature for the day (in ºC)
l.temp.time, h.temp.time – hour of the day when l.temp and h.temp occurred
rain – amount of precipitation (in mm)
ave.wind – average wind speed for the day (in km/h)
gust.wind – maximum wind speed for the day (in km/h)
gust.wind.time – hour of the day when gust.wind occurred
dir.wind – dominant wind direction for the day

Now, let’s move on to Part 2 of this tutorial, where we will start by inspecting the data and prepare it, so we can then proceed to perform EDA.