参考:佐藤信『統計的官能検査法』第19章 (日科技連)
シェッフェの一対比較法(浦の変法)は、比較順序を考慮するもの。
参考文献通りに計算し、(p)値を求める部分だけ追加した。変数名も原則として参考文献のままなので、詳細はそちらを参照してもらいたい。
## -*- coding: utf-8 -*-
## 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"))
p <- pf((S.alpha/df.alpha)/(S.e/df.e), df.alpha, df.e, lower.tail=FALSE)
cat(sprintf("Main | %9.4f %4d %9.4f %9.4f %9.4f", S.alpha, df.alpha, S.alpha/df.alpha, (S.alpha/df.alpha)/(S.e/df.e), p))
if (0.05 < p && p <= 0.1) {
cat(sprintf(" .\n"))
} else if (0.01 < p && p <= 0.05) {
cat(sprintf(" *\n"))
} else if (0.001 < p && p <= 0.01) {
cat(sprintf(" **\n"))
} else if (p <= 0.001) {
cat(sprintf(" ***\n"))
} else {
cat(sprintf("\n"))
}
p <- pf((S.alphaB/df.alphaB)/(S.e/df.e), df.alphaB, df.e, lower.tail=FALSE)
cat(sprintf("Main x Indiv | %9.4f %4d %9.4f %9.4f %9.4f", S.alphaB, df.alphaB, S.alphaB/df.alphaB, (S.alphaB/df.alphaB)/(S.e/df.e), p))
if (0.05 < p && p <= 0.1) {
cat(sprintf(" .\n"))
} else if (0.01 < p && p <= 0.05) {
cat(sprintf(" *\n"))
} else if (0.001 < p && p <= 0.01) {
cat(sprintf(" **\n"))
} else if (p <= 0.001) {
cat(sprintf(" ***\n"))
} else {
cat(sprintf("\n"))
}
p <- pf((S.gamma/df.gamma)/(S.e/df.e), df.gamma, df.e, lower.tail=FALSE)
cat(sprintf("Combi | %9.4f %4d %9.4f %9.4f %9.4f", S.gamma, df.gamma, S.gamma/df.gamma, (S.gamma/df.gamma)/(S.e/df.e), p))
if (0.05 < p && p <= 0.1) {
cat(sprintf(" .\n"))
} else if (0.01 < p && p <= 0.05) {
cat(sprintf(" *\n"))
} else if (0.001 < p && p <= 0.01) {
cat(sprintf(" **\n"))
} else if (p <= 0.001) {
cat(sprintf(" ***\n"))
} else {
cat(sprintf("\n"))
}
p <- pf((S.delta/df.delta)/(S.e/df.e), df.delta, df.e, lower.tail=FALSE)
cat(sprintf("Order | %9.4f %4d %9.4f %9.4f %9.4f", S.delta, df.delta, S.delta/df.delta, (S.delta/df.delta)/(S.e/df.e), p))
if (0.05 < p && p <= 0.1) {
cat(sprintf(" .\n"))
} else if (0.01 < p && p <= 0.05) {
cat(sprintf(" *\n"))
} else if (0.001 < p && p <= 0.01) {
cat(sprintf(" **\n"))
} else if (p <= 0.001) {
cat(sprintf(" ***\n"))
} else {
cat(sprintf("\n"))
}
p <- pf((S.deltaB/df.deltaB)/(S.e/df.e), df.deltaB, df.e, lower.tail=FALSE)
cat(sprintf("Order x Indiv | %9.4f %4d %9.4f %9.4f %9.4f", S.deltaB, df.deltaB, S.deltaB/df.deltaB, (S.deltaB/df.deltaB)/(S.e/df.e), p))
if (0.05 < p && p <= 0.1) {
cat(sprintf(" .\n"))
} else if (0.01 < p && p <= 0.05) {
cat(sprintf(" *\n"))
} else if (0.001 < p && p <= 0.01) {
cat(sprintf(" **\n"))
} else if (p <= 0.001) {
cat(sprintf(" ***\n"))
} else {
cat(sprintf("\n"))
}
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"))
}
3つの試料A/B/Cの場合には、A-B、A-C、B-A、B-C、C-A、C-Bの6通りの一対比較が考えられる。3人による7段階評定((-3)〜(+3))を行ったときには、以下のようにデータを準備して計算する。すると、分散分析表と信頼区間の表が出力される。((p)値より、主効果に有意差があることがわかる)
## usage example (from the reference book)
numStim <- 3
numSubj <- 3
Y <- matrix(data=c(
## comparison result of pairs A-B, A-C, B-A, B-C, C-A, and C-B
-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
## ------------------+-----------+-----------+-----------+----------