close all force
% Y -- Yields of turnips
% pho--Level of phosphate  0 1 2 4 8 16
% lim--Level of liming 0 1 2
%
a=6;   b=3;
Y=[706  1028   928 1185  1120  1050  ...
      998  1025 1111 1111    980  1143 ...
      589    998    961 1014   877  1053];
pho = [ 0 1 2 4 8 16    0 1 2 4 8 16    0 1 2 4 8 16]; 
lim = [ 0 0 0 0 0 0     1 1 1 1 1 1     2 2 2 2 2 2];
%
[p tab stats]=anovan(Y, {pho,lim},'varnames',{'phosphate','liming'},'model','interaction');
Yi=reshape(Y,[a b]);  %Data as a matrix
alphahat = stats.coeffs(2:a+1);  
betahat=stats.coeffs(a+2:a+b+1);
sse=tab{4,2} %  9.3575e+04, SSE from ANOVA table output
d=0; s=0;
for i=1:a
    for j=1:b
        d=d+Yi(i,j)*alphahat(i)*betahat(j);
        s=s+alphahat(i)^2 * betahat(j)^2;
    end
end
% or avoiding the double loop
% d=sum(sum(Yi .* kron( alphahat, betahat' )))
% s =sum(sum(kron( alphahat.^2, betahat'.^2)))
ssab=d^2/s                             %3.2671e+04
ssrem = sse - ssab                  %6.0904e+04
f = ssab/(ssrem/(a*b-a-b))    %4.8280
p = 1-fcdf(f, 1, a*b -a -b)       %0.0556
