Rを用いたシェッフェの一対比較法(中屋の変法)

統計

ワインにうるさい人でも、個別試飲毎に揺るぎない評価ができるかと言えば難しいと思われる。我々、画像を扱う者にとっても撮影条件を変化した画像の良し悪しを評価することは経験上、非常に困難である。しかし、それが2銘柄ずつ(2つに画像)の相対比較であれば絶対評価より容易になる。

このように、一対比較に対してx段階の間隔尺度に基づく相対評価点をさせ、分散分析に持ち込むことで評価対象に有意な差があるかどうかを検定する方法がScheffeの一対比較法である。

Scheffの一対比較法の原法 、たとえばアンケート調査のような場合に用いる。駅前で「アンケートにご協力ください。所要時聞は3 分です。」と頼んでも余程の暇人の協力しか得られないことが予想される。そこで、全対比較の評価を頼むのではなく、一部の一対比較のみをお願いし,その代わりに評価者数を増やすことで検定の信頼性を上げる方法である。

しかしながら、研究室等で主観評価実験を行う場合は、被験者集めを行うことが困難である。そこで一対比較の全組合せを評価してもらうことが多い 。このような場合向けの変法が提案されている。

Scheffの一対比較法の原法と変法の特徴
各被験者が全部対比較評価
しないする
順序効果の考慮あり 原法浦の変法
順序効果の考慮不要芳賀の変法中屋の変法

今回は中屋の変法についてRのコードを作成した。参考にした文献は長沢 伸也 ,川栄 聡史.Excelでできる統計的官能評価法.日科技連出版社である。

Rコード シェッフェの一対比較法(中屋の変法)

## -*- coding: utf-8 -*-

## Scheffe's Pairwise Comparison
## (Nakaya's Variation, not considering stimulus order effect)
## 
## 2012-09-12 by MARUI Atsushi (marui@ms.geidai.ac.jp)
##
## 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"))
  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"))
  }
  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"))
}

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

今回はMRIで撮像した画像から撮影ポジショニングの違いによる画像評価を行った。

画像評価者:3名

評価を行った画像:4枚(A,B,C,D)

A-B, A-C, A-D, B-C, C-Dの総当たりの組み合わせから1対ごとの画像評価を行った。

画像の評価は[-3,-2,-1,0,1,2,3]の7段階評価とした

-3は非常に悪く、3は非常に良い評価となる。

Rコード データセット

numStim <- 4
numSubj <- 3
Y <- matrix(data=c(
## comparison result of pairs A-B, A-C, A-D, B-C, B-D, C-D
               -2, -1, -1,  0, -2,  0,   # subject 1
               -1, -2, -1,  3, 2, 0,   # subject 2
                0,  -1, -1,  0, 0,  2),   # subject 3
              ncol=numSubj, byrow=FALSE)
scheffe.nakaya(Y, numStim, numSubj)

本調査結果を示す。

表と同様にデータを準備して計算すると、分散分析表と信頼区間の表が出力される。(p値より、主効果に有意差があることがわかる)

Rコード 結果

Show in New WindowClear OutputExpand/Collapse Output

> numStim <- 4
> numSubj <- 3
> Y <- matrix(data=c(
+ ## comparison result of pairs A-B, A-C, A-D, B-C, B-D, C-D
+                -2, -1, -1,  0, -2,  0,   # subject 1
+                -1, -2, -1,  3, 2, 0,   # subject 2
+                 0,  -1, -1,  0, 0,  2),   # subject 3
+               ncol=numSubj, byrow=FALSE)
> scheffe.nakaya(Y, numStim, numSubj)

Average Preferences
---------------
a1 =   -0.8333
a2 =    0.5000
a3 =    0.2500
a4 =    0.0833
---------------

ANOVA Table
--------------+---------------------------------------------
Source        |   SS        df   MS         F         p
--------------+---------------------------------------------
Main          |   12.1667    3    4.0556    3.0417    0.1143
Main x Indiv  |   11.3333    6    1.8889    1.4167    0.3415
Combi         |    3.5000    3    1.1667    0.8750    0.5045
Error         |    8.0000    6    1.3333
Total         |   35.0000   18
--------------+---------------------------------------------
                Y(0.05)=1.6319, Y(0.01)=2.3444

Confidence Interval
------------------+-----------------------+----------------------
                  |         95%CI         |         99%CI        
     ai - aj      +-----------+-----------+-----------+----------
                  |   +1.6319 |   -1.6319 |   +2.3444 |   -2.3444 
------------------+-----------+-----------+-----------+----------
a1-a2 =   -1.3333 |    0.2985 |   -2.9652 |    1.0111 |   -3.6778
a1-a3 =   -1.0833 |    0.5485 |   -2.7152 |    1.2611 |   -3.4278
a1-a4 =   -0.9167 |    0.7152 |   -2.5485 |    1.4278 |   -3.2611
a2-a3 =    0.2500 |    1.8819 |   -1.3819 |    2.5944 |   -2.0944
a2-a4 =    0.4167 |    2.0485 |   -1.2152 |    2.7611 |   -1.9278
a3-a4 =    0.1667 |    1.7985 |   -1.4652 |    2.5111 |   -2.1778
------------------+-----------+-----------+-----------+----------

Rを用いたシェッフェの一対比較法(中屋の変法)についてでした。

ご参考までに

コメント

タイトルとURLをコピーしました