シェッフェの一対比較法 (中屋変法) / Scheffe's ANOVA for Paiwise Comparison (Nakaya Variation)

% Scheffe's ANOVA (Nakaya Variation)
%
% Reference: SATOU Shin "Toukei-teki Kannou Kensa-hou (Statistical Affection Test Methods)"
%            (Nikka Giren Publishing, 1985), pp.263--270.
%
% qtukey.m by Trujillo-Ortiz and Hernandez-Walls is necessary to run this m-file:
% http://www.mathworks.com/matlabcentral/fileexchange/3469
%
% 2009-03-23 MARUI Atsushi (marui@ms.geidai.ac.jp)


%% BEGIN "change here"

data = [
% A&B A&C B&C  % compared stimulus-pair (you can increase the number of pairs)
  -2  -1   1   % observer 1
  -2  -2   2   % observer 2
  -3   0   3   % observer 3
  -3   0   2   % observer 4
  -1   1   2   % observer 5
  -3  -1   0   % observer 6
  ];

confInt = 0.95;   % confidence interval level

%% END "change here"




% constants
numSubj = size(data, 1);                  % number of subjects or observers (denoted as N in the book)
numItem = (1+sqrt(1+8*size(data,2)))/2;   % number of items to compare (denoted as t in the book)
%numLevel = 7;                            % number of response levels (not used)





% make this true if you want to know the order the "data" should be placed
if false
  for ni = 1:numItem-1
    for nj = ni+1:numItem
      disp([ni nj]);
    end
  end
end





%% read data
x_ijk = zeros(numItem, numItem, numSubj);
z = 0;
for ni = 1:numItem-1
  for nj = ni+1:numItem
    z = z+1;
    for nk = 1:size(data,1)
      x_ijk(ni,nj,nk) =  data(nk,z);
      x_ijk(nj,ni,nk) = -data(nk,z);
    end
  end
end





%% population estimates
alpha_i = sum(sum(x_ijk, 3), 2) / (numItem*numSubj);   % average preference
alpha_j = sum(sum(x_ijk, 3), 1) / (numItem*numSubj);
alpha_ik = squeeze(sum(x_ijk, 2)/numItem) - repmat(alpha_i, 1, numSubj);   % individual differences
alpha_i_j = repmat(alpha_i,1,numItem) + repmat(alpha_j,numItem,1);
gamma_ij = sum(x_ijk, 3)/numSubj - (alpha_i_j);   % interaction effect





%% sum of squares
S_alpha = sum(sum(sum(x_ijk,3),2).^2) / (numItem*numSubj);
S_alpha_B = sum(sum(sum(x_ijk,2).^2)) / numItem - S_alpha;
S_gamma = sum(sum(triu(sum(x_ijk, 3).^2))) / numSubj - S_alpha;
S_T = sum(sum(sum(x_ijk.^2)))/2;
S_e = S_T - S_alpha - S_alpha_B - S_gamma;





%% degree of freedom
dof_alpha = numItem-1;
dof_alpha_B = (numItem-1)*(numSubj-1);
dof_gamma = (numItem-1)*(numItem-2)/2;
dof_e = (numItem-1)*(numItem-2)*(numSubj-1)/2;
dof_T = numSubj*numItem*(numItem-1)/2;





%% display result
fid = 1;   % fid = fopen('output.txt', 'w');

fprintf(fid, '\n\n\n');
fprintf(fid, '===================================================\n');
fprintf(fid, '=== RESULT OF SCHEFFEs ANOVA (NAKAYA VARIATION) ===\n');
fprintf(fid, '===       (computed on %s)     ===\n', datestr(now, 'yyyy-mm-dd HH:MM:SS'));
fprintf(fid, '===================================================\n\n');

fprintf(fid, '<< POPULATION ESTIMATES >>\n\n');

fprintf(fid, 'alpha_i (mean preferences)\n');
for ni = 1:numItem
  fprintf(fid, '%9.4f\n', alpha_i(ni));
end
fprintf(fid, '\n');

fprintf(fid, 'alpha_ik (individual differences)\n');
for ni = 1:numItem
  for nk = 1:numSubj
    fprintf(fid, '%9.4f', alpha_ik(ni,nk));
  end
  fprintf(fid, '\n');
end
fprintf(fid, '\n');

fprintf(fid, 'gamma_ij (interaction effect)\n');
for ni = 1:numItem
  for nj = 1:numItem
    fprintf(fid, '%9.4f', gamma_ij(ni,nj));
  end
  fprintf(fid, '\n');
end
fprintf(fid, '\n\n');

fprintf(fid, '<< SUM OF SQUARES >>\n\n');

fprintf(fid, 'S_alpha (main effect)\n');
fprintf(fid, '%9.4f\n', S_alpha);
fprintf(fid, '\n');

fprintf(fid, 'S_alpha_B (main effect x individual)\n');
fprintf(fid, '%9.4f\n', S_alpha_B);
fprintf(fid, '\n');

fprintf(fid, 'S_gamma (interaction effect)\n');
fprintf(fid, '%9.4f\n', S_gamma);
fprintf(fid, '\n');

fprintf(fid, 'S_e (error)\n');
fprintf(fid, '%9.4f\n', S_e);
fprintf(fid, '\n\n');

fprintf(fid, '<< ANOVA >>\n\n');

fprintf(fid, '          SS        df   Var       F0\n');
fprintf(fid, '--------- --------- ---- --------- ---------\n');
fprintf(fid, 'S_alpha   %9.4f %4d %9.4f %9.4f\n', S_alpha, dof_alpha, S_alpha/dof_alpha, S_alpha/dof_alpha / (S_e/dof_e));
fprintf(fid, 'S_alpha_B %9.4f %4d %9.4f %9.4f\n', S_alpha_B, dof_alpha_B, S_alpha_B/dof_alpha_B, S_alpha_B/dof_alpha_B / (S_e/dof_e));
fprintf(fid, 'S_gamma   %9.4f %4d %9.4f\n', S_gamma, dof_gamma, S_gamma/dof_gamma);
fprintf(fid, 'S_error   %9.4f %4d %9.4f\n', S_e, dof_e, S_e/dof_e);
fprintf(fid, '--------- --------- ---- --------- ---------\n');
fprintf(fid, 'S_total   %9.4f %4d\n', S_T, dof_T);
fprintf(fid, '\n');





%% compute and display yard-stick
fprintf(fid, 'F(%d, %d; %.2f) = %.4f\n', dof_alpha, dof_e, 1-confInt, finv(confInt, dof_alpha, dof_e));
fprintf(fid, 'F(%d, %d; %.2f) = %.4f\n', dof_alpha_B, dof_e, 1-confInt, finv(confInt, dof_alpha_B, dof_e));
Y = qtukey(dof_e, numItem, confInt) * sqrt((S_e/dof_e) / (numSubj * numItem));
fprintf(fid, 'Y(%.2f) = %.4f\n', 1-confInt, Y);
fprintf(fid, '\n\n');

fprintf(fid, '<< CONFIDENCE INTERVAL >>\n\n');

z = 0;
fprintf(fid, '                   %d%% Confidence Interval\n', round(confInt*100));
fprintf(fid, 'i  j                %+9.4f %+9.4f\n', -Y, Y);
fprintf(fid, '-- --               --------- ---------\n');
for ni = 1:numItem-1
  for nj = ni+1:numItem
    z = z+1;
    fprintf(fid, '%2d %2d =%9.4f  : %9.4f %9.4f\n', ni, nj, alpha_i(ni)-alpha_i(nj), (alpha_i(ni)-alpha_i(nj))-Y, (alpha_i(ni)-alpha_i(nj))+Y);
  end
end
fprintf(fid, '\n\n');




% close file if openned
if fid~=1
  fclose(fid);
end