Lecture 11: Correlation and Linear Regression¶

A gentle introduction to correlation and regression for historians¶

... in Python and Pandas¶

Data Science for Historians (with Python)¶

Created by Kaspar Beelen and Luke Blaxill¶

For the German Historical Institute, London¶

This notebook turns to more advanced and scientifically interesting methods. We will discuss linear regression and later on generalised linear models.

So far, we focused on quantifying and interpreting differences between subgroups in our data (place and gender), more specifically estimating the significance between the means.

Now we turn the analysing the relationship between variables, the extent to which variables are associated with each other. This allows us to model historical phenomena and make predictions.

As in the previous notebook, we highly recommend you watch the lecture by Luke Blaxill for an introduction to the concepts and terminology used in this lecture.

Understanding relations between variables: Lifespan and Wealth in late Victorian London¶

We continue with the data used previously but change the research question and analyse the relationship between wealth and lifespan. Concretely, we will asess if the wealthier districts causes a higher ratio of older people. We don't have information on the possession or income for each individual but we have access to the rateable value per capita for each borough in London, which serves as a proxy for the average prosperity of place. We will abbreviate this variable as rvc in what follows.

When it comes to lifespan we use the same synthetic census data as in notebook ten. We don't have information on how long people live, but we can establish the ratio of 'old' people, for example, those aged over 50.

Notice how both variables (wealth and lifespan) are proxies for concepts we can not measure directly. In many ways, this makes this case study highly problematic (but therefore morer useful from a didactic perspective). However, working with imperfect proxies is often the only option and being creative with your data is a necessity. In the social sciences, the operationalisation of your research is called construct validity: as we project meaning onto our data we should be critical of our semantics.

Let's turn to our code. The first step is to load the required tools we need later on in this notebook.

Before turning to the actual regression analysis, we need to discuss the data wrangling. We will cover a few steps that you will often encounter when preparing your data for analysis.

We create a new dataframe that records the ratio of people over 50 for each district.

Why? We don't know how long Londoners live on average. Here we reason that districts where people have a higher probability of getting "old" (beyond 50) they tend to live longer. For "old", we computed the ratio of residents over 50 years (twice the median) but we could have chosen other thresholds. We encourage you to do this and assess the extent to which the results will change.

To understand the relationship between average wealth and lifespan we need to add additional information to our dataframe.

The dataframe below capture the rateable value per capita in each London borough (these data are synthetic but follow the pattern of "real" values collected for 1922).

We want to merge the data frames, but first, we should check to what extent this is possible, both practically and theoretically.

Our first question is: What variables occur in both datasets? Which ones can we use to combine information from different sources?

For our data, both data frames have information about place. The census lists registration districts while the rateable value is organised by the borough.

This leads to a second question: To what extent are the districts that appear in the census the same as the boroughs in the data on rateable value per capita?

First we want check which place names appear in both datasets. We create a Python set, a list of unique values, for the place names in both dataframes.

We then compute the intersection of these two sets, i.e. the names that appear in both dataframes.

The syntax is set1.intersection(set2). We basically apply the .intersection() to a set and pass another set as argument.

We can merge these two data frame based on the place names. However, be careful: it is not because the strings match that they mean the same thing!

Don't worry about the syntax we will explain it in more detail below.

.merge() is an extremely useful method combining data frames. It requires us to specify which variable (or column) is shared across data frames and can be used for merging tables.

In this .merge() operation, we combined information from rvc with age_by_dist.

As we have shown before, both data frames contain a column with similar values, namely, place names in London.

We use these values to combine the data frames, i.e. rows that have the value Bethnal Green in the age_by_dist data frames will be joined by the information we have for 'Bethnal Green' in rvc

Please notice: the names of the columns can be different but the values in the columns have to overlap.

Below we show a toy example slight adapted from very thorough Pandas documentation.

Imagine that I recorded the number of cats and dogs I spotted each day, but kept this information in two separate data frames.

Both tables have a column that records that day of the week, even though the column name is slightly differentz. We can combine these data frames into one using left_df.merge(right_df). The method is applied to a dataframe at the left-hand side (before the dot) and one at the right-hand side (passed as an argument to the method).

To merge properly, we have to instruct Pandas precisely which columns we want to use for combining information. left_on (right_on) indicates the column name we use in the dataframe at the left-hand (right-hand) side of the .merge() method.

Sometimes we want .merge() using the index of a data frame instead of a column. In this case, we select the argument left_index (or right_index) and which takes a boolean value (True or False).

The .merge() operation should be more understandable at this point: we merge age_by_dist (left) with rvc. For the left data frame we use the index (left_index=True), for the right the 'borough' column (right_on='borough').

Correlation¶

After combining information from these two data frames, we can explore the relationship between variables. Let's start with the most obvious: correlation.

Correlation tells us the extent to which two variables are related. Imagine, each day I collect the number of dogs and cats I observed in two separate variables or vectors (dogs and cats).

If there is a correlation, I expect that if the number of (spotted) dogs up, the same will apply to cats.

Visually, we can make a scatterplot to assess whether there is a correlation: the x-axis shows the number of dogs, the y-axis the number of cats.

The figure suggests a strong correlation between the vectors fordogs and cats, i.e. we observe high values for y when x is high (and vice versa of course).

Please notice that computing correlation is based on an implicit alignment between the two vectors. Each position captures the same day. If the order in cats would be different from dogs (for example the former starts on Wednesday and the latter on Sunday) I could, technically, still make a scatterplot and computer correlation but the results would be meaningless (or at least not match the intention).

The visualization is helpful, but we'd like to put a number on this, quantify the strength of the relation, so we can compare if other vectors are stronger correlated. Here the Pearson correlation coefficient will come to our rescue.

Pearson Correlation¶

A widely used measure for computing the relationship between two vectors or vectors is the Pearson correlation coefficient. This coefficient will vary between -1 and +1, meaning a strong negative respectively positive correlation. A score of 0 indicates no correlation.

We will first show you how to compute the Pearson correlation coefficient yourself, purely for didactic purposes, so you know what is going on (the score is not magic, far from!). But as you can probably guess by this point, there exists a method in Pandas that will help you. We'll have a look at that later.

Showing how the calculation is done, demystifies the process and helps you to understand what is going on.

We need a few ingredients to compute the Pearson correlation coefficient. First, we calculate the mean and standard deviation for each vector. We need to normalize by the amount of data we have and record the number of observations n, which is equal to the length of v1 (or v2—remember both are aligned and equal in length!).

Next, we subtract each value a vector with the mean (of that vector of course). We repeat this for v1 and v2.

Then we compute products of the values in sub_mean_v1 and sub_mean_v2 position-wise, i.e we multiply the first value in sub_mean_v1 with first the value sub_mean_v2 etc.

A convenient function for doing this zip(), which, as the name suggests, "zips" two lists based on their positional index. For example if we zip sub_mean_v1 and sub_mean_v1 we get:

Each element in the list diffs_zipped is a tuple. We multiply the elements in each tuple..

... and we sum all the values, which will give is the nominator of Pearson correlation coefficient.

To obtain the coefficient, we divide the nominator by the product of the standard deviations and the number of observations.

Looking closer at these steps you can intuitively grasp what determines the strength of the correlation. If two values (at the same position) differ from the mean in equal terms, this will result in a higher product (remember from high school that the product of two negative numbers is positive!) and, in turn, increases the value of the nominator.

Of course, Pandas provides a method called.corr() to compute correlations between two instances of pd.Series. We first convert each vector to an instance of pd.Series and then calculate the correlation coefficient, which, no surprise, should be exactly the same!

Of course, now we can also compute the relation between age and wealth.

Actually, that's quite a low correlation! This does not bode well for our expectation that richer boroughs are associated with longer lifespans. But we'll turn that to that shortly.

Spearman Correlation¶

Pearson correlation is one option of many to compute the strength of the relationship between two variables. The choice depends on your research priorities and some procedures will be more relevant or valid than others.

Pearson takes into account the absolute counts of the number of cats and dogs. But imagine, that I don't care about how many animals I observe each day what matters is that ranking of the days is correlated. In this case, the ranking of the days has priority. In this scenario, we should calculate the Spearman rank correlation.

To obtain this coefficient, we need to compute the rank of each day, for example, the third item in v1 has rank 4 (i.e. it comes fourth after the items on positions 0, 1 and 4)

Correlation does not mean causation. Oftentimes, correlations are 'spurious', i.e. it is unlikely there exists a realistic mechanism that explains the association between two variables.

There is even a website on this topic!

Be careful when making claims based on correlation. As often, the statistical experiment won't give you an answer to your question. There is no simple way to determine if a correlation is spurious. It depends on your domain expertise to make sense of these results and justify them.

Linear Regression¶

Whereas correlation tells you that two variables, X and Y, are related you can not model their relationship, i.e. you can not predict Y using X. Correlation can only tell you that when X goes up, Y goes up (or not).

Regression modelling attempts to model the mechanistic relationship between two variables: how many cats do I expect to see, given that I observed five dogs?

Modelling the mechanistic relationship enables us to predict outcomes based on independent variables or predictors.

Returning to our example of age and wealth (goodbye doggies, you were helpful), we'd want to predict the lifespan (the response variable approximated by the ratio of residents older than 50) given the rateable value per capita (predictor variable).

Notice that we rely on domain expertise to construct a model that makes sense. We could also model the reverse, assuming that increasing lifespan will push rateable value upwards.

However, we assumed wealth has a positive impact on lifespan and not vice versa. Our causal model, therefore, looks as follows:

Wealth -> lifespan

and not

Wealth <- lifespan

Whether a causal model makes sense is not a statistical question but depends on your domain expertise as a historian (and in this case, there is room for dispute!?). When reading a historical paper that uses statistical methods, it is often more important to critically look at the relation between the narrative and the method than just the regression procedure: What are data supposed to capture and how are they combined into/framed as a causal theory?

Let's discuss these questions with some actual data. We first demonstrate what a regression model does before we continue with a more detailed explanation of how to run and evaluate a regression model.

A linear model tries to predict a value (lifespan/ratio of 50+) based on a linear combination of other variables (for now, just one, the rateable value per capita). We want to predict the former using the latter.

Using seaborn we can make a dispersion plot where each observation is turned into a point on the figure. The x-axis (x=) shows the rateable value, the y-axis (y=) the ratio of residents older than 50.

A linear regression models the relation between x and y by drawing a line best "fits" the data, i.e. is closest to the observations as possible.

From your high school education, you may remember the equation of a line.

$y = a + b*x$

Which can convert this equation to our variables of interest

• $a$ is the intercept, the point where the line crosses to y-axis when x is zero. You may also think as: what ratio do we expect for a place with zero rateable value?

• $b$ is the slope, and determines the strength of the relationship between $y$ and $x$: how much change do we expect in $y$ after one unit change $x$.

We know lifespan and rateable value.

lifespan = intercept + slope * rateable value per capita

We use these data to estimate the values of $a$ (intercept) and $b$ slope. We will discuss later how, but for now, we just plot the outcome of running a linear regression, to show what the model thinks the line should look like after observing the data.

The line visualizes the outcome of linear regression. The fact that a model throws back a result, a line, doesn't necessarily mean that this result makes sense. The best fit can still be pretty bad. In this case, the line is pretty flat, suggesting that $b% is close to zero, but also pretty bad (it doesn't capture much of the variation).

As you observe there are many values are clustered together in the lower range of the x-axis, and only a few appear at the right-hand side of the figure.

We can improve the fit by transforming the values on x-axis. We take the natural logarithm of the rateable value per capita, this will make small differences bigger and bigger ones smaller, it pulls apart the observation at the lower end of the independent variable and pushes those placed at the far right towards the other observations).

But what is the natural logarithm? Let's have a look, what happens if we take the logarithm of 1000

We get a number close to 7 (i.e. 6.907755278982137). What does this mean? The logarithm is inverse of exponentiation: if x = log(y) then y = exp(x), or sticking to our example, if np.log(1000) produces 6.907755 then e^6.90775 will approximately equal 1000 (not exactly as rounded the numbers. e is called Euler's number and is approximately equal to 2.71828182845904523536028747135266249775724709369995 according to Wikipedia. You can check this for yourself

..is equivalent to...

But why would we do this? Why would we effectively rescale our variables? Taking the logarithm partly solves the problem of dispersion: the logarithmic scale makes linear modelling often easier if a variable is linear on the logarithmic scale.

It is difficult to model these gaps with a line...

However, on a logarithmic scale, these values become neatly fit a line.

Even though 'real' data never fits a pattern as neatly, we observe similar changes when inspecting the rateable value per capita variable.

Both lines don't fit the data perfectly, However, the latter one seems to approximate the actual observed values better.

Also, if now model age as a function of the logged rateable value, the line will fit our observation better across all values of the x-axis. Still the result isn't great...

Later we will quantify goodness of fit of the different lines (before and after taking the logarithm). For now, a visual inspection suffices.

(Please notice that logging also improve the correlation between the two variables.)

Taking the logarithm values doesn't harm, but you'll have to be careful as it influences your interpretation of the results and the predictions of the model.

A closer look at the scatterplot and the data show that two (maybe ever three) places have a significantly higher rateable value: Holborn and Westminster, and also Marylebone. This is confirmed by the boxplot below.

How to approach these observations: one easy, albeit problematic, approach is to consider them as "outliers" and remove them from the datasets.

This is a classical example of where interventions need to be motivated or at least explained by domain expertise.

Let's remove these districts to understand the extent to which they affect our model and more importantly our understanding of the effect of wealth on lifespan. As you'll notice, the figure below suggests that after excluding Holborn and Westminster, the relation become positive. These two districts have a large effect on the slope of the regression line.

The orange line shows the impact of removing the so-called 'outliers'. Throwing away information (i.e. observations) is a tricky thing: yes, it becomes easier to model the relationship between the variables linearly. However, this "solution" comes at a cost, we won't be able to accurately say something about the wealthier districts. If we were to believe the orange line, we'd disastrously overestimate the ratio of older residents in districts like Westminster.

However, if you're only interested in predicting the age for a certain type of boroughs, and have good reasons to believe that the other ones away are qualitatively different (and deserve to be treated as separate case studies) than discarding observations is to some extent justifiable (but, of course, this should be made explicit).

In any case, we need to acknowledge that a few places carry a lot of weight in determining the outcome and interpretation of the regression model (they influence the estimate of the slope substantially).

Again, domain expertise has priority over statistics: the latter alone can not tell you what is the right approach.

Hopefully, this visual introduction demonstrated the use of linear regression and some of the problems that you are likely to encounter. Again, we'd like to stress that often the solution to these issues often require you to lean on your domain expertise as a historian. To solutions are often not statistical.

History always should come before statistics.

The regression model attempts to capture as much as possible variation by drawing a straight line as close as to all the observations. However, the line doesn't have to be straight: this was a choice we made to constrain our model. We could have used, what is called, a "quadratic function" (i.e. $y$ = $a$ + $b$$x$ + $c$$x$²) to obtain a curve. Below we plot the quadratic approximation for the complete (blue) as well as the reduced (orange) data.

Linear Regression (a closer look)¶

In the preceding section, we introduced linear regression from a more visual, intuitive perspective. We did not discuss how to fit and read the results. Let's, therefore, have a closer look at the statistical nitty-gritty of regression analysis. Running a linear regression model is rather straightforward with the statsmodels package in Python. There a just a few things to be careful about. First, we define the response variable, in this case 'Age'.

Next, we select the independent variable from the dataframe.

Now there is one small trick to keep in mind. The equation of a simple linear model is:

$y = a + b*x$

Both $y$ and $x$ are data we observed. We don't know $a$ and $b$, these are the parameters we want to learn from the data. Running a regression will return an estimate for $a$ and $b$ that is most consistent with the observed data.

Because we want statsmodels to estimate both the intercept $a$ and the slope $b$ we have to add one more column with only a constant to X. Luckily this is really straightforwards with the add_constant function.

Now we are ready to fit the model. 'Fitting the model' means finding the optimale fit with the data, in our case we estimate the parameters for a line (intercept and slope) that best capture the variation in our data. Please notice that a line may not necessarily provide

We save the fitted model in the results variable.

Applying the .summary() method to the fitted model will print a report.

The most important parts of this (admittedly extensive) report are located at the top right (R-squared) and in the middle (const and rvc_log). Let's start with the latter. These numbers will figure most prominently in the interpretation of the results.

• const 12.9145 means that for a place with zero rateable value, we expect the percentage of those over 50 to be around 12.9. This is the estimate for the intercept, i.e. the value of y (Age) when x (rateable value) = 0.
• rateable_value_pc 0.0097 is an estimate of the slope which is the increase in y for one unit change in x. For example: if the rateable value changes from 1 to 2 we expect an increase of 0.0097% (not a lot in other words).
• What do the other statistics reported tell us about these coefficients? Continuing with the estimate for rateable_value_pc the p-value is around 0.65, substantially higher than the traditional 0.05 threshold. So far there is a weak positive association, it is unlikely to be of any significance.
• The confidence intervals range from -0.032 to 0.052. Notice that zero is included in the confidence interval. A slope of zero implies a horizontal line, meaning no change in y for x (they are largely independent, x does not influence y at all!). The chance that x and y are not related (in any positive or negative way) is well within the range of possibilities.
• p-values and confidence intervals remain rather abstract and difficult to explain, and statisticians tend to doubt their use. In human language, these number tells you how sure the model is about where to draw the line. It is really sure that the intercept is higher than zero, which is not surprising. However, it is less confident about the slope, which can go down (-0.032) as well as up (0.052).

What happens if we fit a new model with logged rateable values?

Using the log-scale changes coefficients quite a lot (especially for the rateable value). To coefficient is higher (we expect 0.5% change in y for a unit change in x). Nonetheless, the p-value suggests that this positive association may still be the result of random chance. The model is now more certain the relation is positive but still includes zero within the 95% confidence interval.

Note that the scale of the x-axis is different now. For each change in x, the logged rateable value, we expect to add half a percentage to y. The steps on x-axis are become increasingly large on the original scale of rateable value.

According to the model, the rate of change in lifespan is the same when moving from 7 to 20 or transitioning from 20 to 54.

As said, the additional statistics for each coefficient convey where the model thinks the line should be.

Similar to notebook 10, We can make these results more intuitive using a "bootstrap" method. We construct the mean and confidence intervals by resampling from the original data.

More precisely, we repeat the following steps, let's say, 1000 times:

• sample observations from the data (50% at each iteration)
• fit the model and record the values for the intercept and the slope
• plot the lines using these recorded values (or collect them and compute the 2.5 and 97.5 percentile range).

This procedure is implemented by the code below (only 100 iteration but feel free to change this).

The values are not the same because we have a small sample and only ran a few iterations. However, it gives a more intuitive way of understanding what the p-values and confidence intervals of the coefficients mean.

The "bowtie" shape is very common when plotting the outcome of a regression model: the model is fairly certain the line has to go through the mean but becomes more uncertain at the extremes because it lacks data.

Instead of bootstrapping, you can compute the confidence intervals directly with statsmodels. In the figure below we plot the predictions (red) together with the confidence intervals (green).

Please notice that coeficients and their confidence intervals are estimates of the line not the actual data. These are, as you can see, well outside the range given by the predicted values and their confidence intervals.

The .summary() method provides other statistics that help you understanding how well the model performs.

Lastly, we can see what happens when we remove the influential observations or outliers. The estimate for the slope is positive, with a p-value well below 0.05. The model is confident that the association is positive, and generally, knows better where to position the line. This is also reflected in the bootstrap: the lines are tighter packed than in the previous experiments and tend not to go in unpredictable directions.

Before we close this chapter, we need to discuss a more important statistic reported in the regression summary: the R-squared. This score varies between 0 and 1 and indicates how well our model fits the data, with 1 meaning a perfect fit. You'll notice that R-squared in the summary of results_2 is lower compared to the output of results_3 (after removing the outliers).

In general, the R-squared is helpful when comparing models, but we should not attach too much importance to this score: a good fit doesn't necessarily mean that we capture and understand the (causal) mechamism of interest accurately. Firstly, we could simple add more predictors and after while we can almost perfectly fit the data. For example instead of quadratic, a polynomial of order 2, we can fit our data with polynomial of order 20. This will fit the data perfectly but doesn't learn us anything about the relation between our variables of interests. Also the weird squigly lines don't make sense and will make very poor and blatantly unrealistic predictions (a negative percentage -200 for place with a logged rateable value of 2.5. This models overfits the data. Secondly, there is a distinction between prediction and inference. Being good at perdiction, doesn't imply that the model captures the causal mechanism (or the data generation process) adequately. More intelligents thoughts about the difference between prediction and inference are available in Statistical Rethinking or consult this blog post.

This detour aside, we close this notebook with computing R-squared by hand to give you an intuition of the meaning of this statistic.

RMSE: First we compute the root mean squared prediction error. For each data point,

R-squared:

• we compute the variation between observed and predicted values (var_pred)
• we compute the total variation by subtracting the observed value with the mean.

Then we divided var_pred by var_total.

To make sure the number ranges between 0 (minimum) and 1 (maximum), we simply calculate 1 - var.

Fin.¶