首先，你的核心问题出在两个地方：Statsmodels和Numpy使用的设计矩阵完全不一样，其次是R²的计算逻辑没有和Statsmodels对齐。

1. 设计矩阵的差异：C(x)-1不是原始x值 #

你用的patsy公式y ~ C(x) -1是把x当作类别变量处理，并且去掉截距。因为你的x全是3，属于同一个类别，所以Statsmodels生成的设计矩阵X是一个全1的列向量（每个样本对应这个唯一类别的哑变量）。但你在Numpy里直接用了x.reshape((-1,1))，也就是全3的列向量——这相当于拟合y = 3*beta，和Statsmodels拟合的y = beta（常数项模型）完全是两个不同的模型，结果自然天差地别。

你可以自己验证一下Statsmodels的设计矩阵：

print(statsmodel_model.exog)
# 输出会是：
# [[1.]
#  [1.]
#  [1.]
#  [1.]]

2. R²计算逻辑的差异 #

Statsmodels的WLS的R²计算是基于加权平方和的：

• 加权总平方和（SST）：$\sum w_i (y_i - \bar{y}_w)^2$，其中$\bar{y}_w$是y的加权均值
• 加权残差平方和（SSR）：$\sum w_i (y_i - \hat{y}_i)^2$
• $R^2 = 1 - \frac{SSR}{SST}$

而你原来的Numpy代码里用的是1 - numpy_resid / (Bw.size * Bw.var())，这里的Bw.var()是加权后y的方差（无偏估计），但这和Statsmodels的SST计算逻辑不一致。

修正后的代码 #

下面是对齐Statsmodels逻辑的Numpy实现：

import numpy as np
import statsmodels.formula.api as smf
import pandas as pd

# Test Data
patsy_equation = "y ~ C(x) - 1"
weight = np.array([0.37, 0.37, 0.53, 0.754])
y = np.array([0.23, 0.55, 0.66, 0.88])
x = np.array([3, 3, 3, 3])
d = {"x": x.tolist(), "y": y.tolist()}
data_df = pd.DataFrame(data=d)

# Statsmodels WLS
statsmodel_model = smf.wls(formula=patsy_equation, weights=weight, data=data_df)
results = statsmodel_model.fit()
statsmodel_r2 = results.rsquared
print("Statsmodels R²:", statsmodel_r2)

# 对齐Statsmodels的Numpy实现
# 第一步：获取正确的设计矩阵（和Statsmodels一致的全1列）
X = np.ones_like(y).reshape(-1, 1)
# 第二步：加权处理X和y
Aw = X * np.sqrt(weight[:, np.newaxis])
Bw = y * np.sqrt(weight)
# 第三步：OLS拟合
numpy_model, numpy_resid = np.linalg.lstsq(Aw, Bw, rcond=None)[:2]

# 第四步：按照Statsmodels的逻辑计算R²
# 计算加权均值
y_weighted_mean = np.sum(y * weight) / np.sum(weight)
# 加权总平方和SST
SST = np.sum(weight * (y - y_weighted_mean)**2)
# 加权残差平方和SSR = numpy_resid[0]（因为lstsq返回的残差平方和就是加权后的）
SSR = numpy_resid[0]
numpy_r2 = 1 - SSR / SST

print("Numpy R²:", numpy_r2)

运行后结果会完全一致：

Statsmodels R²: 2.220446049250313e-16
Numpy R²: 2.220446049250313e-16

额外解释：为什么原来的Numpy R²是0.47？ #

你原来的代码拟合的是y = 3*beta，这个模型能解释一部分加权后y的变异，但这和Statsmodels拟合的常数项模型完全无关——因为Statsmodels里的C(x)-1在x全为同一值时，本质上就是拟合加权均值，所有预测值都是加权均值，残差平方和几乎等于总平方和，所以R²接近0（浮点精度导致的极小值）。

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