Can Regression Recover the Truth?

Adapted from NUS DSA3361 Tutorial 4.

Suppose we have 100 observations, \((x_1,y_1),\ldots,(x_{100},y_{100})\), and they look roughly like \(y\approx\beta_0+\beta_1x\). You may already know how to use simple linear regression to estimate \(\hat\beta_0\) and \(\hat\beta_1\), giving us the fitted equation

We can use this equation to predict \(y\) for a new \(x\). But have you ever wondered:

Are these computed numbers reliable or not?

Let’s make this question concrete with a familiar example: the formula for a falling object.

Suppose \(y_i\) is the distance an object has fallen, and \(x_i=t_i^2\), where \(t_i\) is the time since it started falling. If you remember a little high-school physics, you may recognise that, ignoring air resistance and starting from rest,

This relationship does not come from statistics or machine learning. It comes from physics. So the true intercept is \(0\), and the true slope is \(\frac12g\), which is about \(4.9\).

But suppose you have forgotten the value of \(g\), and you want to estimate it. We can always grab a small ball, let it fall, and record its falling distance and time. This gives us a dataset. Then we fit an SLR model and obtain the fitted line \(\hat y=\hat\beta_0+\hat\beta_1x\). Since the true slope is \(\frac12g\), twice the fitted slope gives us an estimate of gravity: \(\hat g=2\hat\beta_1\).

Let’s see what this looks like in practice. We will build the dataset step by step: collect times, measure the ideal heights, add inevitable measurement error, and then fit a line.

Try the experiment a few more times and keep track of what you get. Did you notice that \(\hat\beta_1\) changes a little each time?

In short, measurement error in each experiment means that our estimates \(\hat\beta_0\) and \(\hat\beta_1\) will not be exactly equal to the true values \(\beta_0=0\) and \(\beta_1=g/2\). But we still hope that a useful scientific method behaves sensibly: when the measurement error is not too large, OLS—or linear regression more generally—should tend to recover the underlying parameters reasonably well.

That is what we will investigate next. Let’s create a world where we know the truth and see how well regression can recover it.

Generate a world where we know the truth

We generate 100 sample points from

Here’s what each part means:

• \(X_i\) is the observed predictor, generated from \(N(0,1)\).
• \(\epsilon_i\) is a random error term generated from \(N(0,0.25)\).
• The true intercept is \(\beta_0=-1\), and the true slope is \(\beta_1=0.5\).

Let’s create these vectors in Python:

import numpy as np

np.random.seed(123)
x = np.random.normal(loc=0, scale=1, size=100)
eps = np.random.normal(loc=0, scale=np.sqrt(0.25), size=100)
y = -1 + 0.5 * x + eps

With np.random.seed(123), the realised values (rounded to two decimal places) are below:

Applying the true relation \(Y_i=-1+0.5X_i+\epsilon_i\), we obtain the response vector \(y=(y_1,\ldots,y_{100})^\mathsf{T}\). This is the world where we know the truth.

Use the true model to fit the data.

Then we can fit the linear model and print the two estimates:

import statsmodels.api as sm

X = sm.add_constant(x)
model = sm.OLS(y, X).fit()
beta0_hat, beta1_hat = model.params
print(f"beta0_hat = {beta0_hat:.4f}")
print(f"beta1_hat = {beta1_hat:.4f}")

beta0_hat = -1.0095
beta1_hat = 0.4917

So, using the data \((x_i,y_i)_{i=1}^{100}\), our fitted line is \(\hat y=-1.0095+0.4917x\). Looking at the coefficients, it is pretty close to the true relationship \(Y=-1+0.5X+\epsilon\). Nice! Let’s plot the data together with the true relationship and our fitted line.

import matplotlib.pyplot as plt

plt.scatter(x, y, color="black", alpha=0.75, label="Observed data")
y_pred = model.fittedvalues # Compute predicted y
order = np.argsort(x)
x_sorted, y_pred_sorted = x[order], y_pred[order]
plt.plot(x_sorted, y_pred_sorted, color="red", linestyle="--", label="Fitted line")

plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

We can use the function model.summary() below to check the model’s \(R^2\). It also prints quantities with deeper statistical meaning, such as the \(t\)- and \(F\)-test statistics, AIC, and BIC. If you are curious about what all these numbers mean, you are very welcome to take ST3131 — by the end of that course, you will know this whole table inside out 😊.

What if our used model is wrong?

In practice, for a dataset we have collected, who knows the true relationship leh? We usually try several plausible models and compare how they perform. For our example, the true relationship is \(Y=-1+0.5X+\epsilon\). Now suppose we accidentally choose a quadratic model,

and let’s see what happens. The data \(x\) and \(y\) are already available from above, so we can create the extra predictor \(x^2\), fit the model, and inspect its summary:

# Recall: x and y are already generated above.
x1 = x
x2 = x**2
X = sm.add_constant(np.column_stack([x1, x2]))
model2 = sm.OLS(y, X).fit()