Now that we have cleared the dataset and preliminary validated its quality we are ready to start with the exploration. We start with some simple questions.
How many records do we have? How many variables?
What are the variable names? Are they meaningful?
What type is each variable; e.g., numeric, categorical, logical?
How many unique values does each variable have?
What value occurs most frequently, and how often does it occur?
As we discussed earlier, when we are working with spatiotemporal data, there are two types of data. The records that describe the process (runoff_day) and the information about these records (runoff_stations).
After our initial data wrangling is quite easy to answer these questions for runoff_stations. For example, str function can provide us a quick overview of the number of variables, their names and types can be seen with table.
str(runoff_stations)
table(runoff_stations$country) # For one variable
apply(X = runoff_stations, MARGIN = 2, FUN = table) # For all variables
In the case of runoff_day, and in timeseries in general, we need to take a different approach.
Summary Statistics
Although plotting is a necessary first step, in most cases visual inspection offers only limited information about our data. Each time series has 30 to 75 thousand values. We need to somehow derive the information encapsulated in our data to meaningful information. The most straightforward way to do this is by estimating some summary statistics.
By estimating the mean, the standard deviation and the range, we now know how much runoff to expect1 Mean or average is the expected value of the time series. Standard deviation is a measure of spread around the expected value., and how much we shouldn’t2 To be confident about the range, though, we need a really big amount of values.. What we did is to reduce the amount of information from each time series (thousands of values) to four values; mean, standard deviation, minimum and maximum. We must be careful, though. We want to compress and not remove the available information.3 You might have noticed that we have created a new data.table for this information, instead of adding it in runoff_stations. As the volume of analyses grows larger, then we might end up with some huge data objects with “summary information”. We try to avoid this, by keeping relatively small data objects that can be easily combined. To effectively achieve this there should be a common keyword or id between different data objects. In our case, this is sname. Thus, we will keep all information regarding the statistical properties of our dataset in runoff_stats.
For example, looking at the large differences between mean and maximum runoff, one might wonder whether these values are also outliers. A quick way to answer this is by creating a histogram.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
We can see that all histograms4 More about histograms here. have more or less the same shape. The runoff time series appear to have a skewed distribution. It is now obvious that the larger values follows a pattern and, therefore, we should not remove the maximum values as outliers. They are part of the process, even though they are so far away from the mean.5 The statistical measure to estimate this shape feature of the distribution is skewness and can be estimated with skewness function in moments package.
Data Transformation
Descriptive statistics are useful to summarize information, but most of the times this is not enough. Continuous variables, such as river runoff, can have different values, but at the same time be similar. Therefore, we can try to discretize them into meaningful categories.
In our project, we have already estimated runoff means. How can we divide them in three classes; low, medium and high? The obvious answer is to set some criteria, defined as threshold values. Usually, these are relative and subjective, but sometimes quantiles can help.
runoff_stats[, quantile(mean_day)]
## 0% 25% 50% 75% 100%
## 117 1031 1276 2039 2251
Here, we can use the 25% and 75% quantiles to define low and high values respectively.
Now, as we look at our data, it is easier to detect some patterns and see their internal structure. For instance, we observe that there are three subsets in our time series. There are stations in low altitude that belong to large catchments and have high runoff, there are the mid ones, which are the majority and there are the ones that are on or near to the mountains, that belong to smaller catchments and have low runoff. This sounds reasonable, so we can move on.
Discretization is one of the available procedures for transforming our data. Another one, quite necessary in our project is aggregation. Our time series are in daily time scale, but we are interested in changes over monthly, seasonal and annual scales. To do this, we can either sum or average over the coarser scales.6 The choice of function, depends on the process. Summing temperature has no physical meaning. However, summing precipitation has.
In our case, we will use the mean.
runoff_day[, year := year(date)]
runoff_day[, month := month(date)]
runoff_month <- runoff_day[, .(value = mean(value)), by = .(month, year, sname)]
runoff_month[, date := as.Date(paste0(year, '-', month, '-1'))]
By aggregating in other time scales new information can become available. For instance, in runoff data there is seasonality. This is due to the climate cycles in precipitation and temperature (evaporation). We can use boxplots to investigate this.7 More about boxplots here.
We can see that our time series are becoming quite more readable in this scale. We can also observe a drop in runoff in most of the stations after 2000. Could this be the verification of our EDA hypothesis? It is still too early to tell, but we can check whether this shift is common in all time series using a different transformation. We will normalize them by subtracting the mean and dividing by the standard deviation.8 Alternatively, we may use the scale function, which does exactly the same thing.
By rescaling, we can easily compare the fluctuations of all time series and detect any similarities. As expected, most locations are strongly correlated.
Where Are We?
After all these transformations, we are now ready to start exploring our data. However, before doing so, it is a good moment to make a brief stop and check the spatial properties of the data. A quick way to make this happen is by drawing a map.
We observe that the station locations are well-distributed across the river. We can also get a clear idea of their order as water moves from Alps to Baltic. We add this information both to runoff_summary and runoff records.
Since we made this short stop, it is a good opportunity to see another thing. Looking for outliers and missing values is the first line of analyses regarding data robustness. It is a necessary step taken regardless of data type. A more case-specific approach is checking for consistency. This approach relates to the process(es) represented by our data.
In our project, we have river runoff records and some information about the altitude of the measurements and the catchment area they correspond to. We know that larger catchments are expected to have higher runoff.9 This is the information coming from the process.
Now, let’s check if the runoff is consistently higher as catchment size increases.
It looks good. We used a scatter plot to determine whether there is any correlation between runoff and catchment area. The results match our prior belief about how the system should function. This makes us more confident about the future steps of our analysis.
We can save our new datasets and move forward. Who knows, perhaps land is close.
Summary statistics, distributions and histograms: Chapter 2 in Think Stats, by Alen B. Downey.
Assignments
Navigator’s Tasks
Transform runoff_stats to tidy format and then use it to plot mean, median, minimum and maximum for each location as a scaterplot (x axis station, y axis runoff) with different colors and point types for each statistic.
Estimate the skewness and coefficient of variation for each record (a) as a new column in runoff_stats and (b) as a new data.table.
Can you plot each boxplot (facet) of monthly runoff with different fill colour according to the runoff class?
Use boxplot to plot daily runoff per station. What do you observe regarding outliers? Why do you think this happens?
Create your own classes for area and altitude and plot them in a scater plot similar to this one.
Explorer’s Questions
Which is the difference between the median and the 0.5 quantile?
Why the median and the mean are not the same in Rhine runoff?
Do you notice something strange regarding the location of the stations LOBI and REES? Can you think of a possible explanation?
Which were the months, seasons, years with the highest/lowest runoff at each location? Try to present them in comprehensive way. Feel free to improvise!
(Optional) Which is the average distance between each station in km? Which are the two closest and farest adjacent stations?
