読者です 読者をやめる 読者になる 読者になる

シェッフェの一対比較法(浦の変法)

自分の仕事でよく使うシェッフェの一対比較法(浦の変法)を、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"))
}