如何在 MATLAB 中使用 anova1 计算方差分析后的预测响应变量？

Question

我想在 MATLAB 中复制一个使用 SAS 的简单演示（此处概述）。 SAS 代码是

ods rtf file='residual plot SAS output.rtf';
proc print;
run;
proc reg;
model hddt=bloc entry/r p;
output out=new r=residual p=predicted;
run;
proc plot;
plot residual*predicted;
run;
ods rtf close;

（数据位于帖子底部和链接的 pdf 中）

我的 MATLAB 翻译（忽略输出到文件）：

data = readtable('data.txt', "VariableNamesLine" ,1);
[pw,wtbl,stats] = anova1(data.HDDT,data.BLOC);

我有两个问题，一个是小问题，另一个与标题相关：

(1) 该分析看似简单的单向方差分析，但 SAS 和 MATLAB 的确切输出略有不同。

SAS（来自演示指南）：

Source     SS      df     MS       F     Prob>F
-----------------------------------------------
Groups     1.20463    2   0.6023   0.24   0.7893
Error    220.89537   87   2.53903                
Total    222.10000     89                          

MATLAB：

Source     SS      df     MS       F     Prob>F
-----------------------------------------------
Groups     1.067    2   0.53333   0.21   0.8111
Error    221.033   87   2.54061                
Total    222.1     89                          

这可能是由于 SAS 和 MATLAB 中 F 检验期间数据处理的默认值不同（我无法访问 SAS，但过去曾遇到过类似问题的讨论）。我想尽可能复制 SAS 输出。欢迎任何有关如何操作的提示。

(2) SAS 演示根据生成的模型计算响应变量 (hddt) 的预测值。我不知道如何在 MATLAB 中执行此操作，它似乎没有提供遵循方差分析的类似工具。我能找到的最好的就是“统计”工具箱中的“预测”功能，例如

ypred = predict(mdl,Xnew);

但这是最好的方法吗？如果是的话，用什么作为输入？

请注意，SAS 示例中计算的预测变量的值对我来说没有意义，但我可能会在 stats.stackexchange 上问这是什么意思：

感谢您的帮助。

数据（data.txt）：

PLOT BLOC ENTRY HDDT HT LODG YIELD MOIST
3501 1 8 32 84.5 3 46.8 21.5
3502 1 30 32 82.5 3 78.9 19.8
3503 1 7 31 73.5 3 68.5 20.7
3504 1 26 36 64.5 1 64.5 16.7
3505 1 13 33 75.5 2 69.8 16.4
3506 1 18 31 76.0 2 74.5 16.7
3507 1 27 31 78.0 2 85.5 18.0
3508 1 4 34 81.5 4 67.3 16.4
3509 1 29 30 72.0 2 77.4 16.2
3510 1 19 30 80.5 2 55.5 15.3
3511 1 23 34 70.0 4 68.4 15.0
3512 1 5 33 85.0 5 64.2 13.8
3513 1 2 32 87.0 4 71.2 14.7
3514 1 22 35 79.0 4 62.0 13.3
3515 1 25 36 89.5 4 66.0 14.8
3516 1 6 31 83.0 2 77.3 13.6
3517 1 17 33 83.5 3 99.4 16.6
3518 1 3 32 84.0 3 79.9 16.2
3519 1 24 33 82.5 1 83.1 18.0
3520 1 21 32 80.5 1 84.1 19.4
3521 1 20 31 80.5 2 81.5 19.6
3522 1 16 31 75.5 1 65.1 16.0
3523 1 1 32 88.0 3 67.3 19.4
3524 1 14 33 82.0 2 72.9 18.5
3525 1 28 31 79.5 1 89.9 20.6
3526 1 10 31 85.0 3 82.6 20.8
3527 1 11 31 89.0 2 70.0 16.1
3528 1 15 33 75.5 1 77.3 18.9
3529 1 9 33 79.0 2 82.5 16.1
3530 1 12 32 105.0 2 57.4 16.2
3531 2 23 34 72.0 5 69.7 16.5
3532 2 8 32 88.5 3 79.7 16.4
3533 2 28 31 81.0 1 87.6 16.0
3534 2 6 32 81.5 2 70.3 15.8
3535 2 18 30 81.0 2 74.6 13.7
3536 2 2 33 84.5 3 70.7 15.3
3537 2 20 31 86.5 5 79.5 16.6
3538 2 9 32 79.5 3 80.0 13.4
3539 2 13 33 79.0 2 68.7 14.6
3540 2 27 31 73.5 3 78.2 15.4
3541 2 30 32 82.5 2 74.7 16.4
3542 2 21 33 77.0 2 72.5 16.1
3543 2 25 35 84.0 2 72.6 16.0
3544 2 12 32 93.5 3 60.1 16.0
3545 2 16 31 69.5 3 65.2 15.8
3546 2 29 30 72.5 1 82.7 18.4
3547 2 17 33 82.5 4 79.1 17.1
3548 2 14 33 80.0 2 74.5 19.4
3549 2 10 32 76.5 2 81.4 17.2
3550 2 26 36 66.5 1 50.7 22.6
3551 2 5 33 80.0 2 52.9 22.8
3552 2 7 32 77.0 2 75.1 21.3
3553 2 15 33 76.0 1 67.2 19.5
3554 2 24 33 86.5 4 77.9 18.3
3555 2 19 30 84.0 3 72.7 15.9
3556 2 3 33 85.0 3 73.9 17.4
3557 2 1 33 87.0 4 70.6 18.3
3558 2 4 34 86.0 5 70.1 16.0
3559 2 22 35 81.5 6 57.4 16.2
3560 2 11 31 92.0 3 63.3 15.1
3561 3 21 32 83.0 2 83.0 15.6
3562 3 22 34 85.0 6 56.6 13.7
3563 3 14 33 85.0 4 83.9 14.7
3564 3 7 32 86.5 4 80.6 15.5
3565 3 27 31 82.0 2 96.5 14.8
3566 3 6 32 82.5 2 70.4 15.6
3567 3 26 37 73.0 1 72.4 16.2
3568 3 4 35 85.0 4 80.1 16.0
3569 3 25 37 88.5 3 77.6 16.0
3570 3 9 33 88.0 2 80.5 15.6
3571 3 2 32 83.5 2 79.5 16.6
3572 3 16 31 73.5 1 70.0 16.1
3573 3 24 33 83.5 3 85.7 17.8
3574 3 18 31 82.5 2 73.9 17.4
3575 3 10 32 85.5 3 95.7 17.7
3576 3 15 33 78.5 2 78.8 19.9
3577 3 23 35 76.0 5 74.4 19.5
3578 3 17 33 82.0 3 77.3 18.9
3579 3 13 33 77.5 2 72.4 16.2
3580 3 1 33 92.0 4 67.4 16.2
3581 3 30 32 86.0 4 86.5 19.4
3582 3 28 32 75.5 2 89.4 16.7
3583 3 29 30 76.5 1 72.2 19.3
3584 3 11 31 86.5 2 55.6 15.2
3585 3 5 33 80.5 3 63.0 15.5
3586 3 8 32 79.5 3 74.2 14.1
3587 3 12 32 103.0 2 60.9 14.9
3588 3 3 33 83.0 3 72.9 15.7
3589 3 19 30 83.5 3 76.8 14.1
3590 3 20 30 89.0 4 80.4 15.7

Answer 1

anova

而不是

anova1

。对

anova1

的调用相当于 aov = anova(data.HDDT,data.BLOC)

甚至

hddt = data.HDDT;
blocks = data.BLOC;
y = [hddt(blocks==1) hddt(blocks==2) hddt(blocks==3)];
aov = anova(y)

以下也适用

aov = anova(data(:,{'BLOC','HDDT'}),'HDDT')

就像这样

aov = anova(data(:,{'BLOC','HDDT'}),'HDDT~BLOC');

以上所有内容都返回一个方便的结构，其中包含比 

anova1

更多的分析信息，包括残差，可以使用

 检查

aov.Residuals

并且等于（例如对于因子的第一级）

hddt(blocks==1)-mean(hddt(blocks==1))

方差分析项可以“手动”计算如下：

SqDiffG = size(y,1)*(mean(y)-mean(y(:))).^2;
SS_Groups = sum(SqDiffG(:))
SqDiffE = (y-ones(size(y,1),1)*mean(y)).^2
SSE = sum(SqDiffE(:))
SqDiffT = (y-mean(y(:))).^2;
SS_Total = sum(sum(SqDiffT));

与 SAS 结果的差异可能是因为对 SAS 的调用导致计算一组随机值而不是原始输入响应。