0
$\begingroup$

Per Cosma Shalizi, here, equation 11.2, I understood the chi-squared test for regression residuals to involve calculating:

$$\sum_{i=1}^{n} \frac{(y_i - f(x_i, \theta))^2}{\sigma_i^2}$$

which will then be $\chi^2$ distributed with $n - k$ degrees of freedom.

I wanted to test this by simulation but what I get doesn't seem to make sense (I was expecting to reject the null hypothesis in line with the confidence level).

Note on the simulation, performed n_sims times:

  1. Generate a dataset of size n using the true parameter values.
  2. Fit a linear regression and calculate the residuals.
  3. From the residuals attempt to calculate the $\chi^2$ statistic (note I've tried several variations on this to no avail so clearly have some misunderstanding).

Code:

    import matplotlib.pyplot as plt
    import numpy as np    
    import statsmodels.api as sm    
    from scipy import stats
    
    sigma = 1  # true sigma
    beta = 0.5  # true beta
    n = 25  # number of data points
    
    def generate_fit_single_regression(beta, sigma, n):
        X = np.random.uniform(50, 80, size=n)
        y = X * beta + np.random.normal(0, sigma, size=n)
        X_fit = sm.add_constant(X)
        return X, y, sm.OLS(y, X_fit).fit()
    
    df = n - 2
    threshold_95 = stats.chi2(df).ppf(0.975)
    
    n_sims = 10_000
    statistics = []
    for i in range(n_sims):
        np.random.seed(i)
        X, y, results = generate_fit_single_regression(beta, sigma, n)
        predictions = results.predict()
        residuals = y - predictions
        std_residuals = np.std(residuals)
        chi_squared_statistic = np.sum(residuals**2 / np.var(predictions))
        statistics.append(chi_squared_statistic)
    
    plt.figure(figsize=(10, 8))
    plt.hist(statistics, bins=50);
    plt.axvline(threshold_95, color="red", label="$\\chi^2$ threshold");
    plt.legend();
    plt.title("Distribution of computed $\\chi^2$ statistic");

Any help appreciated!

Arnie

enter image description here

$\endgroup$
5
  • $\begingroup$ You seem to confuse $\sigma_i$ with some kind of estimate of it. That's not what Shalizi means. Second, Shalizi explicitly says the degrees of freedom is $1,$ but you have used n - 2 instead: that is, your "$\chi^2$ threshold" is wrong. $\endgroup$ – whuber Apr 8 at 14:33
  • $\begingroup$ Thanks for the reply, @whuber. Re. $\sigma_i$, what should I use instead? Re. degrees of freedom I thought, given we are summing over $n$ $\chi^2(1)$ random variables we would have a $\chi^2(n-k)$ random variable? Clearly I'm probably just highlighting my confusion again at this point so any instruction welcome. $\endgroup$ – Arnold Davidson Apr 8 at 14:46
  • $\begingroup$ Why do you use np.var(predictions) for $\sigma^2$? $\endgroup$ – stefgehrig Apr 8 at 15:45
  • $\begingroup$ Your $\sigma$ seems to be 1 $\endgroup$ – stefgehrig Apr 8 at 15:53
  • $\begingroup$ Perhaps my question should be: what is the correct way to carry out this test? (I'm aware I'm not doing something right.) $\endgroup$ – Arnold Davidson Apr 8 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.