目的

ブラッドリー・テリーのモデル(Bladley-Terry's model)の解を求める。
準対称性の検定を行う（htest クラスのオブジェクトを返す。結果は print.htest による）。
BradleyTerry パッケージに，これに関する BTm という関数がある。
同じ関数名になってしまっていたので，BTM に改名する。

使用法

BTM(x, constant=1, max.rotation=500, epsilon=1e-10, xlab="Score", main="Bradley-Terry model (Paired Comparison)", file="")

引数

x             正方行列（対角成分は何でも良い）
constant      解の制約条件（解の合計値）
max.rotation  収束計算を行う上限回数
epsilon       収束判定値

obj           BTM 関数が返すオブジェクト。クラス名は "htest" と "BTM"
digits        表示桁数

xlab          軸の名称（デフォルトでは Score）
main          グラフのタイトル（デフォルトでは Bradley-Terry model (Paired Comparison)）
file          結果の画像を出力するファイル名（デフォルトでは空文字列で，ファイル出力しない）

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/BTM.R", encoding="euc-jp")

# ブラッドリー・テリーのモデル
# 東京大学教養学部統計学教室編「基礎統計学III　自然科学の統計学」東京大学出版会
# によったが，他のところにあるのとは答えが違う
BTM <- function(x,                                           # 一対比較の結果を表す正方行列
                constant=1,                                     # 解の制約条件（解の合計値）
                max.rotation=500,                               # 収束計算を行う上限回数
                epsilon=1e-10)                                  # 収束判定値
{
        data.name <- deparse(substitute(x))
        method <- "ブラッドリー・テリーのモデル"
        nc <- ncol(x)                                                # 項目数
        stopifnot(nc == nrow(x))                                # 正方行列でないと分析中止
        if (is.null(dimnames(x))) {                             # 項目名がないときは補完する
                labels <- LETTERS[1:nc]
        }
        else if(is.null(colnames(x))) {
                labels <- rownames(x)
        }
        else {
                labels <- colnames(x)
        }
        diag(x) <- 0
        r <- x+t(x)
        yi. <- rowSums(x)
        theta <- rep(constant/nc, nc)
        names(theta) <- labels
        for (i in 1:max.rotation) {
                theta <- yi./colSums(r/outer(theta, theta, "+"))
                theta.new <- theta*(constant/sum(theta))
                converge <- all(abs(theta.new-theta) < epsilon)
                theta <- theta.new
                if (converge) {
                        break
                }
        }
        if (!converge) {
                warning("収束しませんでした。max.rotation を大きくして再実行してみてください")
        }
        s <- !diag(nc)
        expected <- (r*theta/outer(theta, theta, "+"))
        w <- expected[s]
        chi2 <- sum((x[s]-w)^2/w)
        df <- choose(nc-1, 2)
        P <- pchisq(chi2, df, lower.tail=FALSE)
        names(chi2) <- "X-squared"
        names(df) <- "df"
        return(structure(list(statistic=chi2, parameter=df, 
                p.value=P, method=method, data.name=data.name, observed=x, 
                expected=expected, theta=theta), class=c("htest", "BTM")))
}
# summary メソッド
summary.BTM <- function(obj,                                         # BTM オブジェクト
                        digits=5)                                       # 表示する桁数
{
        cat("\nスコア\n\n")
        print(round(obj$theta, digits=digits))
        cat("\nソートされたスコア\n\n")
        print(round(sort(obj$theta), digits=digits))
        cat("\n観察値\n\n")
        print(round(obj$observed, digits=digits))
        cat("\n期待値\n\n")
        print(round(obj$expected, digits=digits))
}
# plot メソッド
plot.BTM <- function(        obj,                                            # BTM オブジェクト
                        xlab="Score",                                   # 結果グラフの横軸名
                        main="Bradley-Terry model (Paired Comparison)", # 結果グラフの表題
                        file="")                                        # 結果グラフをファイル出力するときにファイル名
{
        theta <- obj$theta
        if (file != "") pdf(file, width=540/72, height=160/72, onefile=FALSE)
        plot(theta, rep(0, length(theta)), pch=19, xlab=xlab, main=main, xaxt="n",
                xlim=range(pretty(theta)), ylab="", yaxt="n", ylim=c(0,0.2),
                bty="n", xpd=TRUE)
        text(theta, 0.0, names(theta), pos=3)
        axis(1, pos=0)
        if (file != "") dev.off()
}


使用例

x <- matrix(c(
	0, 14, 7, 13, 16, 18,
	12, 0, 15, 8, 16, 17,
	16, 9, 0, 12, 14, 12,
	12, 17, 12, 0, 13, 7,
	10, 10, 11, 12, 0, 11,
	4, 8, 12, 16, 13, 0
	), byrow=TRUE, nc=6)
rownames(x) <- colnames(x) <- c("Fi", "Br", "Li", "Or", "Ha", "Bu")
(a <- BTM(x))  # 付値しなくてもこれと同じ内容を表示する
summary(a)     # summary メソッドで，追加の情報を表示する
plot(a)        # plot メソッドで，図を描く


出力結果例

> (a <- BTM(x))  # 付値しなくてもこれと同じ内容を表示する

	ブラッドリー・テリーのモデル

data:  x 
X-squared = 22.3877, df = 10, p-value = 0.0132

> summary(a)     # summary メソッドで，追加の情報を表示する

スコア

     Fi      Br      Li      Or      Ha      Bu 
0.19994 0.18805 0.17856 0.16498 0.12984 0.13863 

ソートされたスコア

     Ha      Bu      Or      Li      Br      Fi 
0.12984 0.13863 0.16498 0.17856 0.18805 0.19994 

観察値

   Fi Br Li Or Ha Bu
Fi  0 14  7 13 16 18
Br 12  0 15  8 16 17
Li 16  9  0 12 14 12
Or 12 17 12  0 13  7
Ha 10 10 11 12  0 11
Bu  4  8 12 16 13  0

期待値

         Fi       Br       Li       Or       Ha       Bu
Fi  0.00000 13.39835 12.14938 13.69726 15.76300 12.99201
Br 12.60165  0.00000 12.31045 13.31663 15.38015 14.39112
Li 10.85062 11.68955  0.00000 12.47432 14.47458 13.51094
Or 11.30274 11.68337 11.52568  0.00000 13.98980 12.49840
Ha 10.23700 10.61985 10.52542 11.01020  0.00000 11.60754
Bu  9.00799 10.60888 10.48906 10.50160 12.39246  0.00000

> plot(a)        # plot メソッドで，図を描く



 参考文献
　　東京大学教養学部統計学教室編「基礎統計学III　自然科学の統計学」東京大学出版会