研究で必要になったので、佐藤信『統計的官能検査法』(日科技連 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
以前に作った「浦変法」もあるので、いつか公開します。