Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Exploratory data analysis: wine, heavy rainfall, Titanic

Heinrich Heine University Düsseldorf

We now work through the first steps of exploratory data analysis – as much of it as we have learned so far – on a few examples.

Wine

We begin with an old dataset, which can be downloaded from the UC Irvine Machine Learning Repository, among other places:

Wine. (1991). UCI Machine Learning Repository.

According to the documentation the dataset has 178 rows (called instances here, as they often are) and 14 variables. The first variable is categorical, encoded as an int from {1, 2, 3}; the others are continuous real variables (hence float). The licence is CC-BY-4.0, so we may use the dataset freely, modify it and share our results, provided we name the original source:

Riccardo Leardi, riclea@anchem.unige.it

The dataset comes as two files:

  • wine.data

  • wine.names

Here wine.data is the dataset proper, without any metadata, and wine.names holds the metadata in human-readable form.

From wine.names we take the following in particular:

These data are the results of a chemical analysis of wines grown in the same region in Italy but derived from three different cultivars. The analysis determined the quantities of 13 constituents found in each of the three types of wines.

  • The attributes are

    1. Alcohol

    2. Malic acid

    3. Ash

    4. Alcalinity of ash

    5. Magnesium

    6. Total phenols

    7. Flavanoids

    8. Nonflavanoid phenols

    9. Proanthocyanins

    10. Color intensity

    11. Hue

    12. OD280/OD315 of diluted wines

    13. Proline

  • Number of Instances

    • class 1 59

    • class 2 71

    • class 3 48

Now that we have successfully downloaded the data and read through the metadata that come with them, we can ask what we actually want to do with them.

One obvious task is to find out whether the three classes (that is, the three producers of the wine) can be told apart by the wine’s chemical composition. In statistical terms: do the distributions of the 13 variables differ enough from class to class for the individual data points (the instances) to be assigned to a class on that basis?

We can set this up as a game: if we delete the first variable, the class, can it be reconstructed from the rest? In the simple version we do not delete the class yet, but use it to estimate three different distributions for each of the 13 attributes. Then, for every instance and ignoring its class, we judge which distribution it most likely belongs to, and check whether we were right.

In the advanced version we split the data into training and test data, or form the three (or more) classes without using the variable “class” at all. First, though, let us take a look at the data.

Quite apart from the classes, one thing is interesting straight away: how are the data distributed? We ask that now for each variable on its own.

<class 'numpy.ndarray'> (178, 14)
[1.000e+00 1.423e+01 1.710e+00 2.430e+00 1.560e+01 1.270e+02 2.800e+00
 3.060e+00 2.800e-01 2.290e+00 5.640e+00 1.040e+00 3.920e+00 1.065e+03]
{np.int64(1): np.int64(59), np.int64(2): np.int64(71), np.int64(3): np.int64(48)}
[ True  True  True]
min 11.03 max 14.83 range 3.8000000000000007 std 0.8095429145285168 n 178
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>

We have to admit that this distribution does not look particularly normal. You could almost make out a bimodal distribution – a mixture, that is, of two normal distributions with expected values near 12 and 13.5. Rather than go on guessing from the picture, let us note that we shall need a method for testing such guesses rigorously.

In our example we can press the class into service as an auxiliary variable:

<Figure size 640x480 with 1 Axes>

And indeed, clear differences show at once. The ‘blue’ wine (class 1) has more alcohol on average than the others, the ‘orange’ wine (class 2) less. Let us draw in the corresponding normal densities:

<Figure size 640x480 with 1 Axes>

Unfortunately it is just as clear that the overlaps are so large that between 12 and 14.5 per cent alcohol by volume no sensible prediction of the class is possible.

For completeness we also model the mixture distribution:

<Figure size 640x480 with 1 Axes>

Having seen that the first attribute, the alcohol content, depends on the class, it is natural to examine all the other attributes the same way. So we rewrite the code generically and look at all the variables at once:

<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>

The pictures suggest that, besides attribute 1 (alcohol), attributes 2 (malic acid), 12 (OD280/OD315 of diluted wines, the protein content) and 13 (proline, an amino acid) might also serve to keep the classes apart.

<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>

From this we see that malic acid is a rather poor indicator of the class, whereas protein and proline content together look considerably better suited to prediction than, say, the alcohol content. More precisely: protein content separates class 3 from classes 1 and 2, and proline content separates class 1 from classes 2 and 3. Neither separation is really clean, however, as the plot already shows.

Finally, we look at the relationship between individual attributes with a scatter plot:

<Figure size 640x480 with 1 Axes>
<Figure size 640x480 with 1 Axes>

Whereas alcohol content and malic acid show no obvious relationship, the protein-proline plot has a pronounced white patch: a high proline content goes together with a high protein content. It would now be interesting to compute measures of this relationship, and also to search the specialist literature for an explanation of how it comes about.

To see whether the relationship is only apparent, and perhaps really arises from a dependence on the class, we look at the scatter plot once more, with the classes distinguished by colour:

<Figure size 640x480 with 1 Axes>

And indeed it is quite plain, as it was in the histograms before, that these two attributes taken together separate the classes reasonably well.

What we would like, of course, is a coordinate-independent procedure for separating the classes, so that we do not have to pick two variables by hand. Later, principal component analysis (PCA) will give us one.

Heavy rainfall events (CatRaRE)

The German weather service (DWD) makes data on heavy rainfall events available:

Catalogues of radar-based heavy precipitation events, as a basis for analysing heavy rainfall activity and its damage potential in Germany since 2001

np.void((1, b'W3_Eta_v20', b'W3_Eta_1', 1, 200101020350, 200101020750, b'2001-01-02T03:49:59.000', b'2001-01-02T07:50:00.000', 4, 51, 44.6, 51.0, 100.0, 1, 286, 183, -157462, -4575645, 8315037, b'Baden-W\xfcrttemberg', b'Breisgau-Hochschwarzwal', b'Feldberg (Schwarzwald)', 2.6, 51.5, 38.0, b'28', b'4', 5, 3, 0, 19, b'SWAZF', 19, b'SWAZF', 19, 12, 16.5, 20, 29, 19, 24.8, 30, 6, 620, 21, 0.1, 0.5, 64.9, 0, 0.4, 85, 312, 312, 1187, 585, 1102.3, 1462, -26.1, -311.9, -8.6, 223.5), dtype=[('OBJECTID', '<i8'), ('CatVersion', 'S10'), ('KEY_FIELD', 'S10'), ('Event_ID', '<i8'), ('Start_Time', '<i8'), ('End_Time', '<i8'), ('Date_START', 'S23'), ('Date_END', 'S23'), ('Duration', '<i8'), ('Area', '<i8'), ('Area_LAEA', '<f8'), ('Area_in_G', '<f8'), ('AShare_in_G', '<f8'), ('RRmax_in_G', '<i8'), ('x_RRmax', '<i8'), ('y_RRmax', '<i8'), ('x_RRmaxPRJ', '<i8'), ('y_RRmaxPRJ', '<i8'), ('AGS_RRmax', '<i8'), ('BDL_RRmax', 'S23'), ('LKS_RRmax', 'S23'), ('GMD_RRmax', 'S23'), ('Eta', '<f8'), ('RRmax', '<f8'), ('RRmean', '<f8'), ('Tmax', 'S12'), ('Tmean', 'S12'), ('SRImax', '<i8'), ('SRImean', '<i8'), ('MPP_flag', '<i8'), ('WL_START', '<i8'), ('WLc_START', 'S6'), ('WL_END', '<i8'), ('WLc_END', 'S6'), ('V3_RRmax', '<i8'), ('V3_MIN', '<i8'), ('V3_AVG', '<f8'), ('V3_MAX', '<i8'), ('V4_RRmax', '<i8'), ('V4_MIN', '<i8'), ('V4_AVG', '<f8'), ('V4_MAX', '<i8'), ('BVLK_RRmax', '<i8'), ('BVLK_SUMME', '<i8'), ('BVLK_D1km2', '<i8'), ('SDL_RRmax', '<f8'), ('SDL_GRAD', '<f8'), ('SDL_MAX', '<f8'), ('VSGL_RRmax', '<i8'), ('VSGL_GRAD', '<f8'), ('VSGL_MAX', '<i8'), ('CLC_RRmax', '<i8'), ('CLC_DOMIN', '<i8'), ('SRTM_RRmax', '<i8'), ('SRTM_MIN', '<i8'), ('SRTM_AVG', '<f8'), ('SRTM_MAX', '<i8'), ('TPI_RRmax', '<f8'), ('TPI_MIN', '<f8'), ('TPI_AVG', '<f8'), ('TPI_MAX', '<f8')])
462.0 25.6
5.051172893136403

Now we have to settle on a research question. One question is always obvious: what distribution do the variables have? We can once again look at the bare frequency, that is, at how the points in time are distributed. We could follow it over the years, or average over the years and analyse how the month influences the frequency. The latter suggests itself: heavy rain is a weather phenomenon and so depends on the seasons, and we can expect a clear relationship.

At the moment the start time of a heavy rainfall event is encoded as an int. We turn that to our advantage: make a string out of the int 200101020350 and extract the month from it:

<Figure size 640x480 with 1 Axes>

So the expected relationship shows clearly: in July and June heavy rainfall events are considerably more frequent than in winter.

Next we can look at the location, say at the variable GMD_RRmax and the place "Düsseldorf":

0
b'D\xfcsseldorf'
7
<Figure size 640x480 with 1 Axes>

We see that in Düsseldorf, June and July have so far been the only months with heavy rainfall events.

Seven events are too few for looking at further correlations, so there is little to be gained from staying with Düsseldorf.

To get an impression of the spatial distribution, we look at a scatter plot of the locations of maximum precipitation in the RADKLIM grid (x_RRmax, y_RRmax):

<Figure size 640x480 with 1 Axes>

You can almost make out Germany – and, above all, that the y axis needs rescaling. So we read up on the scaling of the RADKLIM grid and learn that every grid point corresponds to exactly 1×1 km21 \times 1\,\mathrm{km}^2. The plot therefore merely has to be square. To see a little more of the distribution, we also make every point partly transparent. Large clusters then become visible:

<Figure size 1000x1000 with 1 Axes>

This shows that heavy rainfall events are more frequent in the south than in the north. We might also suspect a relationship with the duration (Duration) of the events:

<Figure size 1000x1000 with 1 Axes>

Instead of the suspected relationship we see something else entirely: the duration is mostly a single-digit number, and longer events are rather rare. We can bring in the spatial relationship as well by distinguishing by colour:

<Figure size 1000x1000 with 2 Axes>

Since there is not much to see here (apart from the rarity of longer events), we look at the analogous plot for the mean precipitation:

<Figure size 1000x1000 with 2 Axes>

Here we begin to suspect that the scale is hard to read because the more extreme events (in light blue) take up so much of it. Excluding those events shows more:

<Figure size 1000x1000 with 2 Axes>

But we get the same effect by adjusting the colormap instead:

<Figure size 1000x1000 with 2 Axes>

That will do on heavy rainfall events for now.

Titanic

Now we get to know matplotlib’s object-oriented interface, using the Titanic dataset. On Kaggle it is part of a ‘competition’.

We download the training dataset as titanic.csv.

Problem: ',' in the field 'Name'
array([(0, 3, b'male', 22, 1, 0), (1, 1, b'female', 38, 1, 0), (1, 3, b'female', 26, 0, 0), (1, 1, b'female', 35, 1, 0)], dtype=[('Survived', '<i8'), ('Pclass', '<i8'), ('Sex', 'S10'), ('Age', '<i8'), ('SibSp', '<i8'), ('Parch', '<i8')])
Total survivors: 342
Total passengers: 891
Survival rate: 0.3838383838383838
array([22, 38, 26, 35, 35, -1, 54, 2, 27, 14])
Wrong mean: 22.85409652076319
Correct mean: 29.847605224963715
29.01418439716312 30.425061425061426
<Figure size 1000x1000 with 1 Axes>
<Figure size 1000x1000 with 1 Axes>

It is hardly surprising that in class 3 (the cheapest, in green) the ages tend to be lower than in the other classes. In passing we also notice that class 3 is more frequent than the others (which should perhaps not surprise us either).

To see whether passenger class and survival are related, a scatter plot is enough:

<Figure size 300x300 with 1 Axes>

The scatter plot shows a correlation: passengers in class 3 more often did not survive, and passengers who survived were more often in class 1. It was only visible, however, because we set alpha = 0.005, and it suffers from ignoring the absolute frequencies of the classes (a base rate fallacy).

The correct thing to do here is to compute the probability that a passenger in class 3 survived:

P(Survived=1∣Pclass=3)=P(Survived=1,Pclass=3)P(Pclass=3)P(Survived = 1 | Pclass = 3) = \dfrac{P(Survived = 1, Pclass = 3)}{P(Pclass = 3)}
P( Survived = 1 | Pclass = 1 ) = 0.6296296296296297
P( Survived = 1 | Pclass = 2 ) = 0.47282608695652173
P( Survived = 1 | Pclass = 3 ) = 0.24236252545824846

We see that a first-class passenger was considerably more likely to survive than to die. For second-class passengers the odds are already worse than 50/50, and fewer than a quarter of third-class passengers survived.

It should not surprise us that these probabilities do not sum to 1: they refer to different populations. For the reverse conditional probability, on the other hand, we do expect a sum of 1:

P(Pclass=3∣Survived=1)=P(Survived=1,Pclass=3)P(Survived=1)P(Pclass = 3 | Survived = 1) = \dfrac{P(Survived = 1, Pclass = 3)}{P(Survived = 1)}
P( Pclass = 1 | Survived = 1 ) = 0.39766081871345027
P( Pclass = 2 | Survived = 1 ) = 0.2543859649122807
P( Pclass = 3 | Survived = 1 ) = 0.347953216374269

This shows that survival status certainly does not tell us which class someone travelled in (because the passengers are not spread evenly over the classes). The two directions P(Survived∣Pclass)P(\text{Survived} \mid \text{Pclass}) and P(Pclass∣Survived)P(\text{Pclass} \mid \text{Survived}) are tied together by the class frequencies; exactly how is the business of Bayes’ theorem in the chapter Covariance and dependence, and the chapter Titanic EDA works the same thing through for the variable Sex.

It is always interesting how representative a sample is, because that tells us how far the relationships we have found might generalise. So we now look at whether the variable Sex is correlated with other variables. Admittedly we do not yet know what we are doing, but let that be motivation for learning later what it was we did:

could not convert string to float: np.bytes_(b'male')
-0.003542250865523883 0.9159102368544852

The second return value, pvalue, is a p-value: the probability, for two variables that really are uncorrelated, of seeing a correlation coefficient at least as far from 0 as ours. Here it is 0.92, so an rr like this would be perfectly ordinary for uncorrelated data: the sample cannot tell the tiny deviation from 0 apart from sampling error. That is all the p-value says – in particular, it does not say how likely it is that sex and age are uncorrelated. What a p-value is and is not is settled in the chapter Statistical tests.

PearsonRResult(statistic=np.float64(0.5433513806577556), pvalue=np.float64(1.4060661308798955e-69))

Here we clearly see r>0.5r > 0.5, and a vanishingly small probability pp of observing such an rr for data that are in fact uncorrelated.

PearsonRResult(statistic=np.float64(-0.1319004906047393), pvalue=np.float64(7.850524283810112e-05))
PearsonRResult(statistic=np.float64(-0.003542250865523883), pvalue=np.float64(0.9159102368544852))
PearsonRResult(statistic=np.float64(0.1146308103894224), pvalue=np.float64(0.0006076214735498901))
PearsonRResult(statistic=np.float64(0.24548896009069918), pvalue=np.float64(1.0739326016045421e-13))

So, naively (following the Scipy Stats documentation), we expect the variable SibSp to be correlated with sex. According to the metadata, SibSp is the number of siblings and spouses travelling with the passenger. Age, too, shows a correlation with SibSp:

PearsonRResult(statistic=np.float64(-0.16434928321340586), pvalue=np.float64(8.120291402857697e-07))
<Figure size 500x500 with 1 Axes>