Tutorial 4 · Can Regression Recover the Truth?

← DSA3361 course contents

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_predictors = np.column_stack([x1, x2])

X = sm.add_constant(X_predictors)

model2 = sm.OLS(y, X).fit()

Therefore, the fitted quadratic model is \( \hat y=-0.9925+0.4929x-0.0134x^2.\) The \(t\)-test suggests that the coefficient of \(x^2\), \(\beta_2\), is not significantly different from \(0\) (Curious? ST3131 covers this😊), and the adjusted \(R^2\) does not improve. In other words, adding a more complex term did not meaningfully improve the model.

y_pred2 = model2.fittedvalues

order = np.argsort(x)

plt.scatter(x, y, color="black", alpha=0.75, label="Observed data")

plt.plot(
    x[order],
    y_pred2[order],
    color="red",
    linestyle="--",
    label="Fitted quadratic curve"
)

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

Takeaway: what did we do today?

Today, we started with a world where the true relationship was known:

We then did three things:

• We generated data from this known relationship.

• We fitted the correct linear model and found that \(\hat\beta_0=-1.01\) and \(\hat\beta_1=0.49\) recovered the underlying relationship reasonably well.

• We then added an unnecessary quadratic term \(x^2\). It did not meaningfully improve the model, reminding us that a more complicated model is not automatically a better one.

So, back to the title:

Can regression recover the truth?

In this simulated world, yes. When the true relationship is linear and the noise is not too large, linear regression can recover the underlying relationship reasonably well from the observed data.

More generally, what we did today is called a simulation study. We first decide how the world works, generate data from that known data-generating mechanism, apply statistical methods, and then check how closely they recover the truth. Because the truth is known, we can directly examine quantities such as

But real-data analysis is different: There, nobody tells us the true data-generating mechanism. We therefore cannot directly check whether “the truth” has been recovered, because quantities such as \(\beta\) are unknown.

Instead, when comparing models on real data, we have to rely on quantities that we can actually observe or estimate, such as out-of-sample \(R^2\), MSE, MAE, and other predictive or diagnostic measures. This is often called empirical evaluation or real-data evaluation.

Both simulation studies and empirical evaluations are common ways to evaluate statistical methods in research:

• Simulation studies: we know the truth and ask whether the method can recover it.

• Empirical evaluations: the truth is unknown, so we compare models using observable evidence.

If you’re a TA teaching this tutorial, here are some companion slides you can use alongside these notes.