自分の仕事でよく使うシェッフェの一対比較法（浦の変法）を、Rの勉強がてら実装してみました。佐藤信『統計的官能検査法』（日科技連出版社）
使用例は以下の通り。
numStim <- 3 numSubj <- 3 Y <- matrix(data=c( ## comparison result of pairs 1-2, 1-3, 2-1, 2-3, 3-1, and 3-2 -2, -1, 3, 1, 0, -2, # subject 1 -2, -2, 1, 2, -1, -2, # subject 2 -3, 0, 3, 3, 1, 0), # subject 3 ncol=numSubj, byrow=FALSE) scheffe.ura(Y, numStim, numSubj)
すると平均嗜好度・分散分析表・信頼区間がこんな感じで出てきます。
Average Preferences
---------------
a1 = -0.9444
a2 = 1.3333
a3 = -0.3889
---------------
ANOVA Table
--------------+---------------------------------------------
Source | SS df MS F p
--------------+---------------------------------------------
Main | 50.7778 2 25.3889 27.2836 0.0003
Main x Indiv | 1.2222 4 0.3056 0.3284 0.8515
Combi | 0.0556 1 0.0556 0.0597 0.8131
Order | 0.0556 1 0.0556 0.0597 0.8131
Order x Indiv | 5.4444 2 2.7222 2.9254 0.1113
Error | 7.4444 8 0.9306
Total | 65.0000 18
--------------+---------------------------------------------
Y(0.05)=0.9188, Y(0.01)=1.2813
Confidence Interval
------------------+-----------------------+----------------------
| 95%CI | 99%CI
ai - aj +-----------+-----------+-----------+----------
| +0.9188 | -0.9188 | +1.2813 | -1.2813
------------------+-----------+-----------+-----------+----------
a1-a2 = -2.2778 | -1.3590 | -3.1966 | -0.9965 | -3.5591
a1-a3 = -0.5556 | 0.3633 | -1.4744 | 0.7258 | -1.8369
a2-a3 = 1.7222 | 2.6410 | 0.8034 | 3.0035 | 0.4409
------------------+-----------+-----------+-----------+----------コードは以下の通り。
この関数は単に計算結果を画面に表示するだけなので、計算結果を次の分析やグラフ作成に使おうと思うと面倒です。aov関数のように構造体のような値の返し方がしたいので、どうすればいいか調べてみないと。あと、cat(sprintf(...))のあたりとか、for文を使っているとか、プログラムも汚いですね……。
## Scheffe's Pairwise Comparison ## (Ura's Variation, considering stimulus order effect) ## ## 2012-08-31 by MARUI Atsushi (marui@ms.geidai.ac.jp) ## ## Reference: Chapter 19, Satoh "Statistical Methods in Sensory Tests" Nikkagiren Publishing (1985) scheffe.ura <- function(Y, numStim, numSubj) { ## data preparation (Stim1 x Stim2 x Subj) X <- array(data=0, dim=c(numStim, numStim, numSubj)) k <- 1 for (r in 1:numStim) { for (c in 1:numStim) { if (r != c) { X[r,c,] <- Y[k,] k <- k+1 } } } ## calculate statistics # average preference Xi.. <- apply(X, 1, sum) X.i. <- apply(X, 2, sum) alphai <- (Xi.. - X.i.) / (2*numStim*numSubj) # individual differences Xi.k <- apply(X, c(1,3), sum) X.ik <- apply(X, c(2,3), sum) alphaik <- (Xi.k - X.ik) / (2*numStim) - alphai # combination effect Xij. <- apply(X, c(1,2), sum) Xji. <- t(apply(X, c(1,2), sum)) gammaij <- (Xij. - Xji.) / (2*numSubj) - (matrix(alphai, nrow=length(alphai), ncol=length(alphai), byrow=FALSE) - matrix(alphai, nrow=length(alphai), ncol=length(alphai), byrow=TRUE)) # average order effect delta <- sum(X) / (numStim*(numStim-1)*numSubj) # individual differences of order effect X..k <- apply(X, 3, sum) deltak <- X..k / (numStim*(numStim-1)) - delta ## calculate sum-of-squares # main effect S_\alpha S.alpha <- sum((Xi.. - X.i.)^2) / (2*numStim*numSubj) df.alpha <- numStim-1 # main x subject S_{\alpha(B)} S.alphaB <- sum((Xi.k - X.ik)^2) / (2*numStim) - S.alpha df.alphaB <- (numStim-1)*(numSubj-1) # combination effect S_\gamma tmp <- (Xij. - Xji.)^2 S.gamma <- sum(tmp[upper.tri(tmp, diag=TRUE)]) / (2*numSubj) - S.alpha df.gamma <- (numStim-1)*(numStim-2)/2 # order effect S_\delta S.delta <- sum(X)^2 / (numSubj*numStim*(numStim-1)) df.delta <- 1 # order x subject S_{\delta(B)} S.deltaB <- sum(X..k^2) / (numStim*(numStim-1)) - S.delta df.deltaB <- numSubj-1 # total S_T S.T <- sum(X^2) df.T <- numStim*(numStim-1)*numSubj # error S_e S.e <- S.T - S.alpha - S.alphaB - S.gamma - S.delta - S.deltaB df.e <- numStim*numStim*numSubj - numStim*numStim/2 - 2*numStim*numSubj + 3/2*numStim - 1 ## average preferences cat(sprintf("\n")) cat(sprintf("Average Preferences\n")) cat(sprintf("---------------\n")) for (r in 1:numStim) { cat(sprintf("a%d = %9.4f\n", r, alphai[r])) } cat(sprintf("---------------\n")) ## ANOVA table cat(sprintf("\n")) cat(sprintf("ANOVA Table\n")) cat(sprintf("--------------+---------------------------------------------\n")) cat(sprintf("Source | SS df MS F p\n")) cat(sprintf("--------------+---------------------------------------------\n")) cat(sprintf("Main | %9.4f %4d %9.4f %9.4f %9.4f\n", S.alpha, df.alpha, S.alpha/df.alpha, (S.alpha/df.alpha)/(S.e/df.e), pf((S.alpha/df.alpha)/(S.e/df.e), df.alpha, df.e, lower.tail=FALSE))) cat(sprintf("Main x Indiv | %9.4f %4d %9.4f %9.4f %9.4f\n", S.alphaB, df.alphaB, S.alphaB/df.alphaB, (S.alphaB/df.alphaB)/(S.e/df.e), pf((S.alphaB/df.alphaB)/(S.e/df.e), df.alphaB, df.e, lower.tail=FALSE))) cat(sprintf("Combi | %9.4f %4d %9.4f %9.4f %9.4f\n", S.gamma, df.gamma, S.gamma/df.gamma, (S.gamma/df.gamma)/(S.e/df.e), pf((S.gamma/df.gamma)/(S.e/df.e), df.gamma, df.e, lower.tail=FALSE))) cat(sprintf("Order | %9.4f %4d %9.4f %9.4f %9.4f\n", S.delta, df.delta, S.delta/df.delta, (S.delta/df.delta)/(S.e/df.e), pf((S.delta/df.delta)/(S.e/df.e), df.delta, df.e, lower.tail=FALSE))) cat(sprintf("Order x Indiv | %9.4f %4d %9.4f %9.4f %9.4f\n", S.deltaB, df.deltaB, S.deltaB/df.deltaB, (S.deltaB/df.deltaB)/(S.e/df.e), pf((S.deltaB/df.deltaB)/(S.e/df.e), df.deltaB, df.e, lower.tail=FALSE))) cat(sprintf("Error | %9.4f %4d %9.4f\n", S.e, df.e, S.e/df.e)) cat(sprintf("Total | %9.4f %4d\n", S.T, df.T)) cat(sprintf("--------------+---------------------------------------------\n")) ## calculate yardsticks Y001 <- qtukey(1-0.01, numStim, df.e) * sqrt(S.e/df.e / (2*numSubj*numStim)) Y005 <- qtukey(1-0.05, numStim, df.e) * sqrt(S.e/df.e / (2*numSubj*numStim)) cat(sprintf(" Y(0.05)=%.4f, Y(0.01)=%.4f\n", Y005, Y001)) cat(sprintf("\n")) ## confidence interval cat(sprintf("Confidence Interval\n")) cat(sprintf("------------------+-----------------------+----------------------\n")) cat(sprintf(" | 95%%CI | 99%%CI \n")) cat(sprintf(" ai - aj +-----------+-----------+-----------+----------\n")) cat(sprintf(" | %+9.4f | %+9.4f | %+9.4f | %+9.4f \n", Y005, -Y005, Y001, -Y001)) cat(sprintf("------------------+-----------+-----------+-----------+----------\n")) for (r in 1:(numStim-1)) { for (c in (r+1):numStim) { z <- alphai[r]-alphai[c] cat(sprintf("a%d-a%d = %9.4f | %9.4f | %9.4f | %9.4f | %9.4f\n", r, c, z, z+Y005, z-Y005, z+Y001, z-Y001)) } } cat(sprintf("------------------+-----------+-----------+-----------+----------\n")) cat(sprintf("\n")) }
