在 MATLAB 中比较两个线性回归模型

Question

我想使用 F 统计量比较两个模型的性能。这是一个可重现的示例和预期的结果：

load carbig
tbl = table(Acceleration,Cylinders,Horsepower,MPG);

% Testing separetly both models
mdl1 = fitlm(tbl,'MPG~1+Acceleration+Cylinders+Horsepower');
mdl2 = fitlm(tbl,'MPG~1+Acceleration');

% Comparing both models using the F-test and p-value
numerator = (mdl2.SSE-mdl1.SSE)/(mdl1.NumCoefficients-mdl2.NumCoefficients);
denominator = mdl1.SSE/mdl1.DFE;
F = numerator/denominator;
p = 1-fcdf(F,mdl1.NumCoefficients-mdl2.NumCoefficients,mdl1.DFE);

我们最终得到F = 298.75 和p = 0，表明 mdl1 明显优于 mdl2，由 F 统计量评估。

有没有在不执行两次 fitlm 并进行所有计算的情况下获得 F 和 p 值？

按照@Glen_b 的建议，我尝试运行coefTest，但是该函数的文档记录很差，结果也不是我所期望的。

[p,F] = coefTest(mdl1); % p = 0, F = 262.508  (this F test mdl1 vs constant mdl)
[p,F] = coefTest(mdl1,[0,0,1,1]); % p = 0, F = 57.662 (not sure what this is testing)
[p,F] = coefTest(mdl1,[1,1,0,0]); % p = 0, F = 486.810 (idem)

我相信我应该使用函数[p,F] = coeffTest(mdl1,H,C) 以不同的零假设 (C) 进行测试。但我真的不知道该怎么做，也没有例子。

Answer 1

不，没有。

Fitlm 适合任意模型。在您的情况下，一个带有截距和一个或三个回归量的回归模型。看起来具有三个回归量的模型可以使用来自具有一个回归量的模型的信息，但这仅在模型存在一些限制并且即使这种重叠信息受到限制时才是正确的。 Fitlm 是一个非常通用的框架，可用于任意模型。因此，在共享信息的同时进行多个回归可能会变得非常复杂并且无法实现。 可以为这两个特定模型自己实现这一点。通常这样的线性回归是使用协方差矩阵来解决的：

 Beta = (X' X) ^-1 X' y

如果 X 是变量作为列的数据，y 是目标变量。在这种情况下，您可以重用协方差矩阵的一部分，您只需要较小回归中的列：加速度的变化。由于添加 2 个新变量会为协方差矩阵添加 8 个值，因此您只需节省 1/9 的时间。此外，最重的部分是反转。因此，时间改进非常非常小。

简而言之，只需做两个独立的回归

Answer 2

这个答案是关于比较两个线性回归模型，其中一个模型是另一个模型的受限版本。

简答：

要对估计的系数向量b 的第 3 和第 4 个元素为零的限制进行 F 检验：

[p, F] = coefTest(mdl1, [0, 0, 1, 0; 0, 0, 0, 1]);

进一步说明：

让b 成为我们的估计向量。 b 的线性限制通常以矩阵形式编写：R*b = r。 b 的第 3 和第 4 个元素为零的限制将被写入：

[0, 0, 1, 0    *    b    = [0
 0, 0, 0, 1]                0];

[0, 0, 1, 0; 0, 0, 0, 1] 矩阵是 coefTest 在文档中所称的H 矩阵。

P = coefTest(M,H), with H a numeric matrix having one column for each
    coefficient, performs an F test that H*B=0, where B represents the
    coefficient vector.

加长版

有时使用这种计量经济学例程，最好自己写出来，这样您就知道真正发生了什么。

删除带有 NaN 的行，因为它们只会增加不相关的复杂性：

tbl_dirty = table(Acceleration,Cylinders,Horsepower,MPG);
tbl = tbl_dirty(~any(ismissing(tbl_dirty),2),:);

做估计​​等等...

n = height(tbl);  % number of observations
y = tbl.MPG;
X = [ones(n, 1), tbl.Acceleration, tbl.Cylinders, tbl.Horsepower];
k = size(X,2);     % number of variables (including constant)

b = X \ y;                 % estimate b with least squares
u = y - X * b;             % calculates residuals 
s2 = u' * u / (n - k);     % estimate variance of error term (assuming homoskedasticity, independent observations)
BCOV = inv(X'*X) * s2;     % get covariance matrix of b assuming homoskedasticity of error term etc...
bse = diag(BCOV).^.5;      % standard errors

R = [0, 0, 1, 0;
     0, 0, 0, 1];

r = [0; 0];          % Testing restriction: R * b = r 

num_restrictions = size(R, 1);
F = (R*b - r)'*inv(R * BCOV * R')*(R*b - r) / num_restrictions;   % F-stat (see Hiyashi for reference)

Fp = 1 - fcdf(F, num_restrictions, n - k);  % F p-val

供参考，可以看p。 Hiyashi 著作《计量经济学》的第 65 章。