ドミナンスデータの交互平均法

西里先生の本，「行動科学のためのデータ解析―情報把握に適した方法の利用」の5章からは順序評定データや一対比較データの双対尺度法の話で，データ行列をドミナンス行列に変換してから分析すると良い，と書いてある。
ここでいう分析は交互平均法のことで，交互平均法のコードは前にも書いたから，特に深く追求してませんでした。

最近たまたま指導学生が一対比較データを取ったので，これでやってごらんと本だけ渡して，交互平均法から何からやり始めたんだけど，ドミナンスデータ担ったところで結果が合わないとの相談が。んなわけあるかい，と計算してみたら，おや俺のコードでも結果が合わないw

ということで慌てて復習，再チャレンジしました。
ポイントはp.84にもあるんだけど，

  • 交互平均法の計算につかう行和・列和は，N \times nのドミナンス表の列和と行和にはそれぞれn(n-1)とN(n-1)を入れる」
  • 交互平均法の計算につかう総和はNn(n-1)となる

ということ。そしてこれはテキストにはないことなんだけど，計算途中のベクトルを中心化する「平均値を引く」という作業が要りません。ドミナンス数の行和・列和は0なので，考える必要がないということみたい。

ということで，交互平均法のコードで平均ベクトルを引くところを取り除き(コメントアウトしてます)，行和(mg.row)，列和(mg.row)，総和(mg.t)を計算し直したドミナンス表交互平均法のコードは以下の通りです。

MRA.dominance <- function(dat){
  nc <- ncol(dat)   # n
  nr <- nrow(dat)   # N
  mg.row <- rep(nc * (nc-1),nrow(dat))
  mg.col <- rep(nr * (nc-1),ncol(dat))
  mg.t   <- nr * nc * (nc-1)
  u <- rnorm(nc)
  tmp <- 0
    for( iter in 1:100){
      #algorithm
      v <- u %*% t(dat)
      v <- v / mg.row
#      av <- (mg.row %*% t(v))/mg.t   #Comment OUT!
#      v <- v - (av * rep(1,nr))      #Comment OUT!
      gy <- max(abs(v))
      v <- v/gy
      v
      gy

      u <- v %*% dat
      u <- u / mg.col
#      av <- (mg.col %*% t(u))/mg.t #Comment OUT!
#      u <- u - (av *rep(1,nc))     #Comment OUT!
      gx <- max(abs(u))
      u <- u/gx
      u
      gx
    }   

    eta2 &lt;- gx * gy  # Eigenvalue
    return(eta2)
}

ためにしテキストに載っているデータを使って計算してみましょう。表5.1の，5人の評定者が8人の受験者を順序付ける(1位が1,8位が8)という例。

ord.data <- matrix(c(3,2,3,1,2,
                     6,7,5,8,7,
                     8,5,7,6,4,
                     1,3,2,3,5,
                     4,1,1,2,1,
                     5,8,6,7,8,
                     2,4,4,5,3,
                     7,6,8,4,6),nrow=5)
ord.data
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    3    6    8    1    4    5    2    7
## [2,]    2    7    5    3    1    8    4    6
## [3,]    3    5    7    2    1    6    4    8
## [4,]    1    8    6    3    2    7    5    4
## [5,]    2    7    4    5    1    8    3    6

これをドミナンス表に変える関数も作ったよ。

dominance <- function(dat){
  nc <- ncol(dat)
  nr <- nrow(dat)
  dom <- matrix(ncol=ncol(dat),nrow=nrow(dat))
  for(i in 1:nc){
    for(j in 1:nr){
      dom[j,i] <- nc+1-2*dat[j,i]
    }
  }
  return(dom)
}

(d.data <- dominance(ord.data))

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    3   -3   -7    7    1   -1    5   -5
## [2,]    5   -5   -1    3    7   -7    1   -3
## [3,]    3   -1   -5    5    7   -3    1   -7
## [4,]    7   -7   -3    3    5   -5   -1    1
## [5,]    5   -5    1   -1    7   -7    3   -3

そしてこれを上で作った関数に入れます。

MRA.dominance(d.data)
## [1] 0.3199945

はい，第一固有値がでました。

という感じです。