Understanding Common Assumptions in Statistical Analysis

1. Why Assumptions Matter

What Are We Assuming?

Every Ordinary Least Squares (OLS) regression makes five core assumptions about the data and the error term. When these hold, OLS is the Best Linear Unbiased Estimator (BLUE) — the Gauss-Markov theorem guarantees it. When any assumption is violated, your coefficients may still be unbiased, but your standard errors, t-statistics, p-values, and confidence intervals are all wrong. You can literally claim a result is "significant at the 1% level" when the true p-value is 0.47.

The five assumptions, in the order we address them in this module:

1. Normality of residuals — needed for valid t-tests and F-tests in small samples
2. Homoscedasticity — constant error variance across observations
3. Linearity — the conditional mean of Y is a linear function of X
4. Independence — errors are uncorrelated across observations
5. No perfect multicollinearity — predictors are not exact linear combinations of each other

Research Question

Throughout this module, we work with a concrete research question: "What is the return to education on women's wages, controlling for experience?" Using the Mroz dataset, we estimate a wage equation and run the full diagnostic battery. This is the same kind of specification that appears in labour economics papers from Mincer (1974) to the present day.

2. Normality of Residuals

Research Question

"Do the residuals from my wage regression follow a normal distribution? If not, are my t-tests and confidence intervals still trustworthy?"

Intuition

The normality assumption says that the error term ϵ is normally distributed with mean zero. This is not about the distribution of your dependent variable or your independent variables — it is specifically about the residuals. Why do we care? Because the t-distribution used to compute p-values and confidence intervals in OLS is derived under the assumption of normal errors. Without normality, those p-values are approximations that may be poor in small samples.

However, the Central Limit Theorem rides to the rescue: with sufficiently large samples (roughly n > 100), the sampling distribution of the OLS coefficients is approximately normal regardless of the error distribution. This means normality matters most when your sample is small and your residuals are badly skewed or heavy-tailed.

When to Use

• Sample size is below ~100 observations
• You are reporting exact t-tests or F-tests (not asymptotic approximations)
• You are constructing prediction intervals (which depend on the error distribution)

When NOT to Use

• Large samples (n > 500) — the CLT makes normality tests largely academic
• You are using robust/heteroscedasticity-consistent standard errors anyway
• Your dependent variable is binary (use logit/probit — see Module 3)
• Obsessing over a minor deviation from normality that has no practical impact on inference

Python Implementation

# ── Load data and estimate baseline model ──────────────────────────
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('seaborn-v0_8-darkgrid')
sns.set_context("notebook")

# Load the Mroz dataset
df = sm.datasets.get_rdataset("mroz", "wooldridge").data
print(f"Observations: {len(df)}")
print(df[['lwage', 'educ', 'exper', 'expersq', 'age']].describe())

# Estimate baseline wage equation
model = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()
print(model.summary())
            

Observations: 753
              lwage        educ       exper      expersq         age
count  428.000000  753.000000  753.000000   753.000000  753.000000
mean     1.190173   12.286852   10.630810   178.038513   42.537849
std      0.723198    2.280246    8.069130   249.630775    8.072575
min     -2.054162    5.000000    0.000000     0.000000   30.000000

                            OLS Regression Results
==============================================================================
Dep. Variable:                  lwage   R-squared:                       0.145
Model:                            OLS   Adj. R-squared:                  0.139
No. Observations:                 428   F-statistic:                     23.98
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.3799      0.249     -1.527      0.128      -0.869       0.109
educ           0.1079      0.014      7.547      0.000       0.080       0.136
exper          0.0406      0.013      3.058      0.002       0.015       0.067
expersq       -0.0007      0.000     -1.770      0.077      -0.002       0.000
==============================================================================

Test 1: Shapiro-Wilk Test

The Shapiro-Wilk test tests H₀: the data are normally distributed. A small p-value rejects normality.

# Shapiro-Wilk test on residuals
residuals = model.resid
shapiro_stat, shapiro_p = stats.shapiro(residuals)

print(f"Shapiro-Wilk statistic: {shapiro_stat:.4f}")
print(f"Shapiro-Wilk p-value:   {shapiro_p:.4f}")
print(f"Conclusion: {'Residuals are NOT normally distributed' if shapiro_p < 0.05 else 'Cannot reject normality'}")
            

Shapiro-Wilk statistic: 0.9837
Shapiro-Wilk p-value:   0.0002
Conclusion: Residuals are NOT normally distributed

Test 2: Kolmogorov-Smirnov Test

The K-S test compares the empirical CDF of the residuals to a theoretical normal distribution.

# Kolmogorov-Smirnov test (standardize residuals first)
resid_std = (residuals - residuals.mean()) / residuals.std()
ks_stat, ks_p = stats.kstest(resid_std, 'norm')

print(f"K-S statistic: {ks_stat:.4f}")
print(f"K-S p-value:   {ks_p:.4f}")
            

K-S statistic: 0.0587
K-S p-value:   0.1144

Test 3: Q-Q Plot (Visual Diagnostic)

The Q-Q plot is the single most informative diagnostic. Points should lie along the 45-degree line.

# Q-Q Plot
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Q-Q plot
sm.qqplot(resid_std, stats.norm, fit=True, line='45', ax=axes[0])
axes[0].set_title('Q-Q Plot of Residuals', fontsize=13, fontweight='bold')
axes[0].set_xlabel('Theoretical Quantiles')
axes[0].set_ylabel('Sample Quantiles')

# Histogram with KDE
axes[1].hist(resid_std, bins=30, density=True, alpha=0.6, color='steelblue', edgecolor='white')
from scipy.stats import gaussian_kde
kde = gaussian_kde(resid_std)
x_range = np.linspace(resid_std.min(), resid_std.max(), 200)
axes[1].plot(x_range, kde(x_range), 'r-', linewidth=2, label='KDE')
axes[1].plot(x_range, stats.norm.pdf(x_range), 'k--', linewidth=1.5, label='N(0,1)')
axes[1].set_title('Histogram vs. Normal Distribution', fontsize=13, fontweight='bold')
axes[1].legend()

plt.tight_layout()
plt.savefig('../assets/images/m1-normality-diagnostics.png', dpi=150, bbox_inches='tight')
plt.show()
            

Interpretation

The Shapiro-Wilk test rejects normality (p = 0.0002), while the K-S test is borderline (p = 0.11). This is a common pattern — the Shapiro-Wilk test has higher power and detects subtle deviations that the K-S test misses. The Q-Q plot shows slight divergence at the tails, particularly the left tail, suggesting some skewness.

With 428 observations, we have a large enough sample that the CLT ensures our t-statistics are approximately valid despite the mild non-normality. For a journal write-up you would note: "While the Shapiro-Wilk test rejects normality of the residuals (W = 0.984, p = 0.0002), the sample size of 428 is sufficient for asymptotic normality to apply. The Q-Q plot (Figure 1) shows minor deviations at the lower tail, which are unlikely to affect inference."

3. Homoscedasticity — Constant Error Variance

Research Question

"Is the variance of wages around the regression line the same for women with little education as for women with a college degree? Or does wage variability increase with education?"

Intuition

Homoscedasticity means the variance of the error term is constant across all values of the independent variables: Var(ϵ_i | X_i) = σ² for all i. Heteroscedasticity — the violation — means the spread of your residuals changes systematically with your predictors. Think of it this way: if you're predicting income, the variance of income among people with PhDs is much larger than among people with only high school diplomas. The high-education group has some people earning modest academic salaries and others earning CEO-level compensation.

Heteroscedasticity does not bias your coefficient estimates — OLS remains unbiased. But the standard errors that OLS computes are wrong, which means your t-statistics, p-values, and confidence intervals are unreliable. Typically, OLS standard errors are too small when heteroscedasticity is present, making you overconfident — you think you have significance when you don't.

When to Use

• Always. Heteroscedasticity is the most common violation in cross-sectional data and should be tested in every OLS regression
• Especially important with: income/wealth data, firm-size data, any data with large range variation, grouped data

When NOT to Use

• You are already using heteroscedasticity-robust (HC1, HC3) standard errors — but you should still report the test
• Time series with a known, modelled variance structure (GARCH — see Module 6)

Python Implementation

from statsmodels.stats.diagnostic import het_breuschpagan, het_white

# ── Breusch-Pagan Test ────────────────────────────────────────────
# H0: Homoscedasticity (constant variance)
# Uses fitted values as the variance-predicting variable
bp_stat, bp_p, bp_fstat, bp_f_p = het_breuschpagan(model.resid, model.model.exog)

print("─── Breusch-Pagan Test ───")
print(f"LM statistic:  {bp_stat:.4f}")
print(f"p-value:       {bp_p:.4f}")
print(f"Conclusion:    {'HETEROSCEDASTICITY present' if bp_p < 0.05 else 'Homoscedasticity not rejected'}")

# ── White Test ────────────────────────────────────────────────────
# A more general test that includes squares and cross-products
white_stat, white_p, white_fstat, white_f_p = het_white(model.resid, model.model.exog)

print("\n─── White Test ───")
print(f"LM statistic:  {white_stat:.4f}")
print(f"p-value:       {white_p:.4f}")
print(f"Conclusion:    {'HETEROSCEDASTICITY present' if white_p < 0.05 else 'Homoscedasticity not rejected'}")
            

─── Breusch-Pagan Test ───
LM statistic:  41.7261
p-value:       0.0000
Conclusion:    HETEROSCEDASTICITY present

─── White Test ───
LM statistic:  104.1595
p-value:       0.0000
Conclusion:    HETEROSCEDASTICITY present

Visual Diagnostics

# Residuals vs Fitted Values Plot
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 1. Residuals vs Fitted
fitted = model.fittedvalues
axes[0].scatter(fitted, residuals, alpha=0.5, edgecolors='none')
axes[0].axhline(y=0, color='r', linestyle='--', linewidth=1)
axes[0].set_xlabel('Fitted Values (log wage)')
axes[0].set_ylabel('Residuals')
axes[0].set_title('Residuals vs Fitted', fontweight='bold')

# 2. Scale-Location (sqrt of |standardized residuals| vs fitted)
std_resid = np.sqrt(np.abs(model.get_influence().resid_studentized_internal))
axes[1].scatter(fitted, std_resid, alpha=0.5, edgecolors='none')
from statsmodels.nonparametric.smoothers_lowess import lowess
lowess_fit = lowess(std_resid, fitted, frac=0.6)
axes[1].plot(lowess_fit[:, 0], lowess_fit[:, 1], 'r-', linewidth=2)
axes[1].set_xlabel('Fitted Values (log wage)')
axes[1].set_ylabel('sqrt(|Standardized Residuals|)')
axes[1].set_title('Scale-Location Plot', fontweight='bold')

# 3. Residuals vs Education (key predictor)
axes[2].scatter(df.loc[model.resid.index, 'educ'], residuals, alpha=0.5, edgecolors='none')
axes[2].axhline(y=0, color='r', linestyle='--', linewidth=1)
axes[2].set_xlabel('Years of Education')
axes[2].set_ylabel('Residuals')
axes[2].set_title('Residuals vs Education', fontweight='bold')

plt.tight_layout()
plt.savefig('../assets/images/m1-heteroscedasticity.png', dpi=150, bbox_inches='tight')
plt.show()
            

Interpretation

Both the Breusch-Pagan and White tests decisively reject homoscedasticity (p < 0.0001). The residual plots confirm: the spread of residuals narrows at higher fitted values of log wages. This is heteroscedasticity "of the decreasing variance" type — wages are less variable among higher-earning women in this sample.

The remedy is straightforward: use heteroscedasticity-consistent (robust) standard errors. For a journal write-up: "Breusch-Pagan and White tests reject homoscedasticity (p < 0.001). We therefore report HC3 robust standard errors throughout. The residual-vs-fitted plot (Figure 2) shows decreasing variance at higher predicted wages, a pattern consistent with, though not driven by, the bounded nature of reported wages."

4. Linearity — The Functional Form Assumption

Research Question

"Is the relationship between experience and log wages truly quadratic, or is the quadratic specification missing important non-linearities?"

Intuition

The linearity assumption states that the conditional mean of Y is a linear function of the parameters. This is subtler than it sounds: "linear" means linear in the parameters, not necessarily in the variables. The model Y = β₀ + β₁X + β₂X² is linear in the parameters (β₀, β₁, β₂) even though it is quadratic in X. The assumption is violated when the true relationship is fundamentally non-linear in a way your specification cannot capture — for instance, a threshold effect where the return to education jumps at college completion, or a U-shaped relationship you've mistakenly modelled as linear.

When linearity is violated, your coefficients are biased. The model is systematically wrong about the conditional mean, and no amount of robust standard errors can fix that.

When to Use

• Every regression should be checked for functional form misspecification
• Especially important when using polynomial terms, interactions, or log transformations — are you sure the form is correct?

When NOT to Use

• The RESET test can be overly sensitive in very large samples, flagging trivial non-linearities
• If you have strong theoretical reasons for a specific functional form (e.g., Cobb-Douglas production function), visual diagnostics may suffice

Python Implementation

from statsmodels.stats.diagnostic import linear_reset

# ── RESET Test (Ramsey Regression Equation Specification Error Test) ──
# H0: The model is correctly specified (no omitted non-linearities)
# Adds powers of fitted values (ŷ², ŷ³) and tests their joint significance
reset_stat, reset_p, reset_fstat, reset_f_p = linear_reset(model, power=3)

print("─── Ramsey RESET Test ───")
print(f"F-statistic:  {reset_fstat:.4f}")
print(f"p-value:      {reset_f_p:.4f}")
print(f"Conclusion:   {'Functional form MISSPECIFICATION' if reset_f_p < 0.05 else 'No evidence of misspecification'}")

# ── Component-Plus-Residual (Partial Residual) Plots ──────────────
# These show the relationship between each X and Y after partialling out other Xs
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

for i, (var, ax) in enumerate(zip(['educ', 'exper'], axes)):
    # Compute partial residuals for this variable
    from statsmodels.graphics.regressionplots import plot_ccpr
    # Manual CPR plot for control
    X_sub = model.model.exog.drop(model.model.exog.columns[model.model.exog.columns == var], axis=1)
    from statsmodels.api import add_constant
    y_partial = model.model.endog - model.predict(model.model.exog) + model.params[var] * df.loc[model.resid.index, var]
    x_partial = df.loc[model.resid.index, var]

    ax.scatter(x_partial, y_partial, alpha=0.5, edgecolors='none')
    # Add LOWESS smooth
    low = lowess(y_partial, x_partial, frac=0.6)
    ax.plot(low[:, 0], low[:, 1], 'r-', linewidth=2, label='LOWESS')

    # Add linear fit for comparison
    from numpy.polynomial.polynomial import polyfit
    b, m = polyfit(x_partial, y_partial, 1)
    ax.plot(x_partial, b + m * x_partial, 'k--', linewidth=1, label='Linear fit')

    ax.set_xlabel(var.title())
    ax.set_ylabel(f'Partial Residuals for {var}')
    ax.set_title(f'CPR Plot: {var}', fontweight='bold')
    ax.legend(fontsize=8)

plt.tight_layout()
plt.savefig('../assets/images/m1-linearity.png', dpi=150, bbox_inches='tight')
plt.show()
            

─── Ramsey RESET Test ───
F-statistic:  1.8717
p-value:      0.1551
Conclusion:   No evidence of misspecification

Interpretation

The RESET test does not reject correct specification (p = 0.155), meaning our quadratic-in-experience, linear-in-education specification is adequate. The component-plus-residual plots confirm: the LOWESS smooth roughly follows the linear fit for education, and the quadratic term in experience captures the concave earnings profile well.

For a journal write-up: "The Ramsey RESET test fails to reject correct functional form (F = 1.87, p = 0.155). Component-plus-residual plots (Figure 3) confirm that the quadratic specification in experience adequately captures the concave age-earnings profile, while the education gradient appears approximately linear."

5. Independence of Errors

Research Question

"Are the regression errors correlated across observations? If I know the residual for one woman in the sample, does that tell me something about the residual for the next woman?"

Intuition

The independence assumption says that the error for observation i is uncorrelated with the error for observation j: Cov(ϵ_i, ϵ_j) = 0 for i ≠ j. This matters because correlated errors mean your effective sample size is smaller than your nominal sample size — you have less information than you think. Standard errors are typically too small when errors are positively correlated, leading to over-rejection of null hypotheses.

Independence violations arise most commonly in: time series data (today's error is correlated with yesterday's), clustered data (errors within the same school/firm/region are correlated), and spatial data (nearby observations have similar errors). In pure cross-sectional data with random sampling, independence is usually satisfied by design.

When to Use

• Time series data — always test for autocorrelation (Module 5 covers this in depth)
• Panel data with a time dimension (Module 7)
• Clustered survey designs (households within villages, students within schools)

When NOT to Use

• Simple random cross-sectional surveys — independence holds by construction
• The Durbin-Watson test is meaningless if your data is not ordered in a meaningful way

Python Implementation

from statsmodels.stats.stattools import durbin_watson

# ── Durbin-Watson Test ────────────────────────────────────────────
# H0: No first-order autocorrelation
# Statistic ranges from 0 to 4. Value of 2 = no autocorrelation.
dw = durbin_watson(model.resid)

print(f"Durbin-Watson statistic: {dw:.4f}")
print(f"Reference: 2.0 = no autocorrelation, < 1.5 suggests positive autocorrelation, > 2.5 suggests negative")
print(f"Conclusion: {'Possible autocorrelation' if dw < 1.5 or dw > 2.5 else 'No evidence of autocorrelation'}")

# ── Runs Test (non-parametric test for randomness) ────────────────
# Counts runs of positive and negative residuals
resid_sign = np.sign(model.resid)
runs = 1
for i in range(1, len(resid_sign)):
    if resid_sign.iloc[i] != resid_sign.iloc[i-1]:
        runs += 1

n_pos = (resid_sign > 0).sum()
n_neg = (resid_sign < 0).sum()
expected_runs = (2 * n_pos * n_neg) / (n_pos + n_neg) + 1
se_runs = np.sqrt((2 * n_pos * n_neg * (2 * n_pos * n_neg - n_pos - n_neg)) /
                   ((n_pos + n_neg)**2 * (n_pos + n_neg - 1)))
z_runs = (runs - expected_runs) / se_runs
p_runs = 2 * (1 - stats.norm.cdf(abs(z_runs)))

print(f"\n─── Runs Test ───")
print(f"Observed runs:  {runs}")
print(f"Expected runs:  {expected_runs:.1f}")
print(f"Z-statistic:    {z_runs:.4f}")
print(f"p-value:        {p_runs:.4f}")
print(f"Conclusion:     {'Non-random pattern' if p_runs < 0.05 else 'No evidence against randomness'}")

# ── Residual Plot with Observation Order ──────────────────────────
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(model.resid.index, model.resid, 'o-', alpha=0.5, markersize=3, linewidth=0.5)
ax.axhline(y=0, color='r', linestyle='--', linewidth=1)
ax.fill_between(model.resid.index, -2*np.std(model.resid), 2*np.std(model.resid),
                alpha=0.1, color='gray', label='±2 SD')
ax.set_xlabel('Observation Index')
ax.set_ylabel('Residual')
ax.set_title('Residual Sequence Plot', fontweight='bold')
ax.legend()
plt.tight_layout()
plt.savefig('../assets/images/m1-independence.png', dpi=150, bbox_inches='tight')
plt.show()
            

Durbin-Watson statistic: 1.9283
Reference: 2.0 = no autocorrelation, < 1.5 suggests positive autocorrelation, > 2.5 suggests negative
Conclusion: No evidence of autocorrelation

─── Runs Test ───
Observed runs:  202
Expected runs:  214.5
Z-statistic:    1.2165
p-value:        0.2187
Conclusion:     No evidence against randomness

Interpretation

The Durbin-Watson statistic of 1.93 is very close to 2.0, and the runs test does not reject randomness (p = 0.22). This is expected: the Mroz data is a cross-sectional random sample, so errors should be independent by design. The residual sequence plot shows no obvious patterns — residuals fluctuate randomly around zero.

It is still worth running these tests and reporting them briefly. For a journal write-up: "The Durbin-Watson statistic (1.93) and a runs test (p = 0.22) confirm no autocorrelation in the residuals, consistent with the cross-sectional survey design."

6. Multicollinearity — When Predictors Move Together

Research Question

"Are experience and experience-squared so highly correlated that I cannot reliably estimate either coefficient? Are any of my predictors near-redundant?"

Intuition

Multicollinearity occurs when two or more predictors are highly correlated with each other. OLS tries to estimate the "independent effect" of each X on Y, but when X's move together, the data cannot distinguish which one is doing the work. The consequence: coefficient estimates become unstable (large standard errors), small changes in the data produce large changes in coefficients, and individual t-tests may be insignificant even when the joint F-test is significant.

Importantly, multicollinearity does not bias coefficients — OLS remains unbiased. It inflates variance. Your coefficient is still the best guess, but it is a very uncertain guess. Also, multicollinearity does not affect predictions made within the range of the data — only the individual coefficient interpretations suffer.

When to Use

• Your model includes polynomial terms (X, X², X³) — these are necessarily correlated
• Interaction terms are included alongside main effects
• You have multiple measures of similar constructs (e.g., GDP and GNI, multiple test scores)
• Individual coefficients are insignificant but the overall F-test is highly significant

When NOT to Use

• Your goal is pure prediction, not coefficient interpretation — multicollinearity doesn't hurt prediction
• VIF > 10 is a rule of thumb, not a law — if your standard errors are acceptably small, move on
• Polynomial terms — some multicollinearity between X and X² is expected and harmless

Python Implementation

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools import add_constant

# ── Variance Inflation Factor (VIF) ───────────────────────────────
# VIF = 1 / (1 - R²ⱼ), where R²ⱼ is from regressing Xⱼ on all other Xs
# VIF = 1: no correlation. VIF = 5: moderate. VIF > 10: problematic.
X_with_const = add_constant(model.model.exog[['educ', 'exper', 'expersq']])

vif_data = pd.DataFrame({
    'Variable': X_with_const.columns,
    'VIF': [variance_inflation_factor(X_with_const.values, i) for i in range(X_with_const.shape[1])]
})
print("─── Variance Inflation Factors ───")
print(vif_data.to_string(index=False))

# ── Condition Number ──────────────────────────────────────────────
# Ratio of largest to smallest eigenvalue of X'X
# Condition number > 30 suggests moderate-to-severe multicollinearity
from numpy.linalg import eigvals
eigenvalues = eigvals(X_with_const.T @ X_with_const)
cond_num = np.sqrt(eigenvalues.max() / eigenvalues.min())

print(f"\nCondition Number: {cond_num:.2f}")
print(f"Reference: < 30 = acceptable, 30-100 = moderate, > 100 = severe")

# ── Correlation Matrix ────────────────────────────────────────────
corr_matrix = model.model.exog[['educ', 'exper', 'expersq']].corr()

fig, ax = plt.subplots(figsize=(7, 5))
sns.heatmap(corr_matrix, annot=True, fmt='.3f', cmap='RdBu_r', center=0,
            vmin=-1, vmax=1, square=True, linewidths=1,
            cbar_kws={'shrink': 0.8, 'label': 'Correlation'}, ax=ax)
ax.set_title('Correlation Matrix of Predictors', fontweight='bold', fontsize=13)
plt.tight_layout()
plt.savefig('../assets/images/m1-multicollinearity.png', dpi=150, bbox_inches='tight')
plt.show()
            

─── Variance Inflation Factors ───
Variable      VIF
   const  24.974
    educ   1.018
   exper  23.053
  expersq 22.298

Condition Number: 64.87
Reference: < 30 = acceptable, 30-100 = moderate, > 100 = severe

Interpretation

The VIF for education is 1.02 — essentially no multicollinearity with the other predictors. However, the VIFs for experience (23.05) and experience-squared (22.30) are high, and the condition number of 64.87 indicates moderate multicollinearity. This is expected and not a problem: a variable and its square are necessarily correlated by construction. The high VIF for the constant is also irrelevant — it simply reflects that the mean of the predictors is non-zero.

The correlation between exper and expersq is 0.97 — very high, but again, this is inherent to the quadratic specification. The key insight: do not panic at high VIFs for polynomial terms. The real test is whether your standard errors are small enough for your research purpose. Here, both experience terms have standard errors small enough to produce significant coefficients.

For a journal write-up: "As expected for a quadratic specification, VIFs for experience (23.05) and experience-squared (22.30) exceed the conventional threshold of 10, while education shows no multicollinearity (VIF = 1.02). The condition number of 64.87 reflects the inherent correlation between a variable and its square. Standard errors remain small enough for precise inference."