Statsmodels部分识别模型回归问题：为何与Stata输出不一致？

Great question! Let's break down what's happening here, and how to get the result you expect.

What's Causing the Unexpected Statsmodels Output? #

The core difference lies in how Statsmodels and Stata handle rank-deficient (perfectly collinear) regression models:

1. Stata's Behavior:
   Stata automatically detects perfectly collinear variables and drops redundant ones to leave a full-rank model. In your case, x1, x2, and the constant term are all identical across both observations (they're always 1), so Stata drops two of them, leaving only x3 (plus maybe the constant, depending on settings). This gives you the expected coefficient of 1 for x3, and a warning that standard errors can't be estimated (since you have exactly as many parameters as observations, leaving 0 degrees of freedom).

2. Statsmodels' Behavior:
   Statsmodels' OLS uses the Moore-Penrose pseudoinverse (generalized inverse) to solve rank-deficient models instead of dropping variables. This method returns the minimum-norm solution—the set of coefficients where the sum of squared coefficients is as small as possible.

   For your data, the model equations are:

   2 = const + x1 + x2 + x3
   1 = const + x1 + x2 + 0

   Subtracting gives x3 = 1, and we're left with const + x1 + x2 = 1. The minimum-norm solution here is const=0, x1=0.5, x2=0.5 (since 0² + 0.5² + 0.5² + 1² = 1.5, which is smaller than any other valid combination like const=1, x1=0, x2=0 whose squared sum is 2). The infinite standard errors make sense too—these coefficients aren't uniquely identified (there are infinitely many valid solutions), so their variance is undefined.

How to Get the Expected Output in Statsmodels #

To replicate Stata's behavior, you need to manually drop redundant collinear variables to create a full-rank model. Here's how:

Step 1: Identify Redundant Variables #

First, check the rank of your design matrix to confirm collinearity:

import numpy as np
import pandas as pd
import statsmodels.api as sm

data = np.array([[2, 1, 1, 1], [1, 1, 1, 0]])
df = pd.DataFrame(data, columns=['y', 'x1', 'x2', 'x3'])
z = df.pop('y')
X = sm.add_constant(df)

# Check rank of the design matrix
print(f"Design matrix rank: {X.rank()}")  # Output: 2
print(f"Number of columns: {X.shape[1]}")  # Output: 4

This tells us we have 2 redundant variables (since rank < number of columns). In your case, const, x1, and x2 are perfectly collinear (all values are 1 across observations).

Step 2: Drop Redundant Variables and Refit #

Choose a subset of variables that form a full-rank matrix. For example, keep only const and x3:

# Keep only constant and x3
X_reduced = X[['const', 'x3']]
mod_reduced = sm.OLS(z, X_reduced)
res_reduced = mod_reduced.fit()

print(res_reduced.summary())

Output Explanation #

The result will show:

• const coefficient = 1.0
• x3 coefficient = 1.0 (matches your expectation)
• Standard errors marked as NaN (since we have 2 observations and 2 parameters, leaving 0 degrees of freedom—just like Stata's warning that standard errors can't be estimated)

Alternatively, you could drop the constant and keep x3 plus one of x1/x2—you'll still get x3=1 with undefined standard errors.

内容的提问来源于stack exchange，提问作者FooBar