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

先日の浦の変法に続き、シェッフェの一対比較法(中屋の変法)を実装してみました。佐藤信『統計的官能検査法』(日科技連出版社)に書かれている方法をそのままコードにして、今度はp値の計算結果も表示するようにしただけです。詳しいことは同書を参照して下さい。

使用例は以下の通り。

numStim <- 3
numSubj <- 6
Y <- matrix(data=c(
## comparison result of pairs 1-2, 1-3, 2-3
              -2, -1,  1,   # subject 1
              -2, -2,  2,   # subject 2
              -3,  0,  3,   # subject 3
              -3,  0,  2,   # subject 4
              -1,  1,  2,   # subject 5
              -3, -1,  0),  # subject 6
              ncol=numSubj, byrow=FALSE)
scheffe.nakaya(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   43.1132    0.0007
Main x Indiv  |   11.2222   10    1.1222    1.9057    0.2470
Combi         |    0.0556    1    0.0556    0.0943    0.7711
Error         |    2.9444    5    0.5889
Total         |   65.0000   18
--------------+---------------------------------------------
                Y(0.05)=0.8323, Y(0.01)=1.2617

Confidence Interval
------------------+-----------------------+----------------------
                  |         95%CI         |         99%CI        
     ai - aj      +-----------+-----------+-----------+----------
                  |   +0.8323 |   -0.8323 |   +1.2617 |   -1.2617 
------------------+-----------+-----------+-----------+----------
a1-a2 =   -2.2778 |   -1.4454 |   -3.1101 |   -1.0160 |   -3.5395
a1-a3 =   -0.5556 |    0.2768 |   -1.3879 |    0.7062 |   -1.8173
a2-a3 =    1.7222 |    2.5546 |    0.8899 |    2.9840 |    0.4605
------------------+-----------+-----------+-----------+----------

コードは以下の通り。

## Scheffe's Pairwise Comparison
## (Nakaya's Variation, not considering stimulus order effect)
##
## Reference: Chapter 20, Satoh "Statistical Methods in Sensory Tests"
## Nikkagiren Publishing (1985)


scheffe.nakaya <- 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-1)) {
    for (c in (r+1):numStim) {
      X[r,c,] <-  Y[k,]
      X[c,r,] <- -Y[k,]
      k <- k+1
    }
  }
  
  ## calculate statistics
                                        # average preference
  Xi.. <- apply(X, 1, sum)
  alphai <- Xi.. / (numStim*numSubj)

                                        # individual differences
  Xi.k <- apply(X, c(1,3), sum)
  alphaik <- Xi.k / numStim - alphai

                                        # combination effect
  Xij. <- apply(X, c(1,2), sum)
  gammaij <- Xij. / numSubj - (matrix(alphai, nrow=length(alphai), ncol=length(alphai), byrow=FALSE) - matrix(alphai, nrow=length(alphai), ncol=length(alphai), byrow=TRUE))
  
  
  ## calculate sum-of-squares
                                        # main effect S_\alpha
  S.alpha <- sum(Xi..^2) / (numStim*numSubj)
  df.alpha <- numStim-1
                                        # main x subject S_{\alpha(B)}
  S.alphaB <- sum(Xi.k^2) / numStim - S.alpha
  df.alphaB <- (numStim-1)*(numSubj-1)
                                        # combination effect S_\gamma
  tmp <- Xij.^2
  S.gamma <- sum(tmp[upper.tri(tmp, diag=TRUE)]) / numSubj - S.alpha
  df.gamma <- (numStim-1)*(numStim-2)/2
                                        # total S_T
  S.T <- sum(X^2)/2
  df.T <- numSubj*numStim*(numStim-1)/2
                                        # error S_e
  S.e <- S.T - S.alpha - S.alphaB - S.gamma
  df.e <- (numStim-1)*(numStim-2)*(numSubj-1)/2

  ## 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("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 / (numSubj*numStim))
  Y005 <- qtukey(1-0.05, numStim, df.e) * sqrt(S.e/df.e / (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"))
}