regression - 在 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 值？

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

[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

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

简答:

要对估计的系数向量 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 页。