2011-11-06

単調増加の検出と比較（８）Isotonic programming回帰結果の評価について

単調増加の検出と比較の目次
前までの記事で、単調減少(増加)へのあてはめが、順序・半順序を問わずできるようになった
さて、これを、多件数に実施して、それらを一塊として評価することを考えよう
IsoGeneパッケージは、「順序」の場合のIsotonic programmingをした上で、その当て嵌まりの良さを評価し、多件数で実行する仕組みである
- 当て嵌まりの良さの評価の部分を踏襲できると思われるので、関数の中をのぞいてみる
- IsoGene1()関数では、次のことをしている
  - 増加パターンか減少パターンかの判定
  - 増加・減少それぞれのパターンについて、実測値が推定値からどれくらいずれているかを定量する
    - 定量方法として５つの方法を用いてすべて算出している
      - Esquare
      - Williams
      - Marcus
      - M
      - ModM
IsoGene1()関数の中身を見る

library(IsoGene)
IsoGene1

とすれば中身が見える(日本語コメントは、この記事のために書き込んだもの)

> IsoGene1
function (x, y) 
{
    ordx <- order(x)
    x <- x[ordx]
    y <- y[ordx]
    unx <- unique(x)
    n.p <- table(x)
    n.g <- length(n.p)
# 群ごとの平均
    y.m <- tapply(y, as.factor(x), mean)
    y.m.tot <- rep(mean(y), length(unx))
    y.is.u <- pava(y = y.m, w = n.p)
    y.is.d <- pava(y = y.m, w = n.p, decreasing = TRUE)
# ここまでが、pava()関数を用いた、『Isotonic modelへの当て嵌め』
# 以降が、実測データと当て嵌め値とのずれを数値化する部分
# u,d, up/dnがup,downを表す
# u/dの向きが決まっていれば、この区別は不要になる
    iso.u <- rep.iso.u <- rep.iso.d <- y.m.all <- NULL
    rep.iso.d <- rep(y.is.d, n.p)
    rep.iso.u <- rep(y.is.u, n.p)
    y.m.all <- rep(y.m, n.p)
# 全レコードの平方和(すべてのレコードで同じ値が返ってくるとしたら…)
    SST0 <- sum((y - mean(y))^2)
# up/downモデルからのずれの平方和
    SSIS.u1 <- sum((rep.iso.u - y)^2)
    SSIS.d1 <- sum((rep.iso.d - y)^2)
# 群ごとの平均とレコードとの作る平方和
    SST <- sum((y - y.m.all)^2)
    direction <- if (SSIS.u1 <= SSIS.d1) 
        "u"
    else "d"
# 以下は５つの方法のそれぞれの流儀で値を算出
    lambda1.up <- SSIS.u1/SST0
    Esquare.up <- 1 - lambda1.up
    iso.u <- y.is.u
    w.up <- (y.is.u[n.g] - y.m[1])/sqrt(sum((y - y.m.all)^2)/(sum(n.p) - 
        n.g) * (1/n.p[1] + 1/n.p[n.g]))
    w.c.up <- (y.is.u[n.g] - y.is.u[1])/sqrt(sum((y - y.m.all)^2)/(sum(n.p) - 
        n.g) * (1/n.p[1] + 1/n.p[n.g]))
    m.up <- (y.is.u[n.g] - y.is.u[1])/sqrt(SSIS.u1/(sum(n.p) - 
        n.g))
    i.up <- (y.is.u[n.g] - y.is.u[1])/sqrt(SSIS.u1/(sum(n.p) - 
        length(unique(y.is.u))))
    lambda1.dn <- SSIS.d1/SST0
    Esquare.dn <- 1 - lambda1.dn
    iso.u <- y.is.d
    w.dn <- (y.is.d[n.g] - y.m[1])/sqrt(sum((y - y.m.all)^2)/(sum(n.p) - 
        n.g) * (1/n.p[1] + 1/n.p[n.g]))
    w.c.dn <- (y.is.d[n.g] - y.is.d[1])/sqrt(sum((y - y.m.all)^2)/(sum(n.p) - 
        n.g) * (1/n.p[1] + 1/n.p[n.g]))
    m.dn <- (y.is.d[n.g] - y.is.d[1])/sqrt(SSIS.d1/(sum(n.p) - 
        n.g))
    i.dn <- (y.is.d[n.g] - y.is.d[1])/sqrt(SSIS.d1/(sum(n.p) - 
        length(unique(y.is.d))))
    res <- list(E2.up = Esquare.up, Williams.up = as.numeric(w.up), 
        Marcus.up = as.numeric(w.c.up), M.up = as.numeric(m.up), 
        ModM.up = as.numeric(i.up), E2.dn = Esquare.dn, Williams.dn = as.numeric(w.dn), 
        Marcus.dn = as.numeric(w.c.dn), M.dn = as.numeric(m.dn), 
        ModM.dn = as.numeric(i.dn), direction = direction)
    return(res)
}
<environment: namespace:IsoGene>