シェッフェの一対比較法(中屋変法)のMatlabコード

研究で必要になったので、佐藤信統計的官能検査法』(日科技連 1985年)を参考に、シェッフェの一対比較法(中屋変法)による分散分析を行うMatlabスクリプトを作りました。実のところ「参考にした」どころではなく、書いてあった手計算のアルゴリズムをそのまま実装したような感じです。

使い方は「BEGIN "change here"」から「END "change here"」の間にあるデータを書き換えて実行するだけ。分散分析表と信頼区間の表を出力します。データ形式は、各行が実験参加者一人一人で、一対比較の評価結果を並べる形にします。一行の中の評価結果の順番は、1&2、1&3、2&4・・・1&t、2&3、2&4・・・2&t、3&4・・・、(t-2)&(t-1)、(t-1)&t、となります。ここで1&2は試料1と試料2の比較結果を表しています。たとえば「菓子Aと比較して菓子Bはどのくらい○○であったか」に-3〜+3の7段階で答えてもらい、全組み合わせの結果を並べて書きます。

% 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


%% BEGIN "change here"

data = [
% 1&2 1&3 2&3
  -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

%% 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

以前に作った「浦変法」もあるので、いつか公開します。