wilcox_test (exact test) の使い方

formula で，~ の右に来る変数は factor でないといけない。よくわからないエラーメッセージでしばらく悩んだ。

> library(coin)
> d <- data.frame(x=c(3,1,2,4,3,2,5,4,3,4,5),
+                 g=c(1,1,1,1,1,2,2,2,2,2,2))
> wilcox_test(x~g, data=d)
 以下にエラー check(itp) : ‘object’ does not represent a two sample problem
> wilcox_test(x~g, data=d, distribution="exact")
 以下にエラー check(itp) : ‘object’ does not represent a two sample problem
> d$g <- factor(d$g)
> wilcox_test(x~g, data=d)

    Asymptotic Wilcoxon Mann-Whitney Rank Sum Test

data:  x by g (1, 2)
Z = -1.5884, p-value = 0.1122
alternative hypothesis: true mu is not equal to 0

> wilcox_test(x~g, data=d, distribution="exact")

    Exact Wilcoxon Mann-Whitney Rank Sum Test

data:  x by g (1, 2)
Z = -1.5884, p-value = 0.1515
alternative hypothesis: true mu is not equal to 0

> wilcox.test(x~g, data=d)

    Wilcoxon rank sum test with continuity correction

data:  x by g
W = 6.5, p-value = 0.1349
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ：
In wilcox.test.default(x = c(3, 1, 2, 4, 3), y = c(2, 5, 4, 3, 4,  :
   タイがあるため、正確な p 値を計算することができません

変数値に対応する要素

変数 x のとる値が，1，3，6，9 のとき，対応する色 red, green, blue, brown を割り当てる
x <- c(1,6,9,3,3,6,9,1,3)
y <- as.integer(names(table(x)))
color <- c("red", "green", "blue", "brown")
sapply(x, function(z) color[which(z==y)]) # 解1
color[matrix(outer(x, y, "=="), length(x)) %*% 1:length(y)] # 解2

プログラムの最適化？

元のプログラム

system.time({set.seed(12345)
hits <- res <- 1:1000
for (j in 1:1000) {
for (i in 1:1000) { res[i] <- ifelse((20 %in% sample(0:99,50,replace=TRUE)),1,0) }
hits[j] <- sum(res) }
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過
#    57.292      0.784     57.598

少しずつ書き換えてみる

system.time({set.seed(12345)
hits <- replicate(1000,
sum(replicate(1000, 20 %in% sample(0:99, 50, replace=TRUE))))
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過
#    30.168      0.651     30.599

次，

system.time({set.seed(12345)
x <- array(sample(0:99, 50*1000*1000, replace=TRUE), dim=c(50, 1000, 1000))
hits <- apply(x, 3, function(y) sum(apply(y, 2, function(z) any(z == 20))))
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過
#    17.678      1.631     19.786

次

system.time({set.seed(12345)
hits <- replicate(1000,  {
x <- matrix(sample(0:99, 50*1000, replace=TRUE), 50)
sum(apply(x, 2, function(y) any(y == 20)))
})
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過
#    18.812      0.450     19.171

次

system.time({set.seed(12345)
x <- matrix(sample(0:99, 50*1000*1000, replace=TRUE), 50)
hits <- colSums(matrix(apply(x, 2, function(y) any(y == 20)), 1000))
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過
#    23.416      1.489     25.095

次

system.time({set.seed(12345)
x <- matrix(sample(0:99, 50*1000*1000, replace=TRUE) == 20, 50)
hits <- colSums(matrix(colSums(x) > 0, 1000))
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過 
#     2.978      3.711     21.218

ユーザ時間は短くなったけど，経過時間は余り差がない（どういうことなんだろうか）

system.time({set.seed(12345)
x <- matrix(sample(0:99, 50*1000*1000, replace=TRUE) == 20, 50)
gc() # ゴミ掃除
hits <- colSums(matrix(colSums(x), 1000) > 0)
hist(hits)
print(summary(hits))
})
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  349.0   384.0   395.0   394.7   405.0   440.0
#   ユーザ   システム       経過 
#     2.818      1.268      4.403

ゴミ掃除やると早くなった

因子間相関係数行列がへん?

いつの間にか（R 2.13.0 ですでに），factanal(   rotation="promax")  で，因子間相関係数行列が出力されるようになっているが，因子との対応がおかしいのではないか?
Loadings は SS loadings の大きい順に並べ替えられているのに，因子間相関係数行列の行と列はその並べ替えができていないのじゃないか?
つまり，
temp <- factanal(x, factors=n, rotation="none")
ans <- promax(loadings(temp))
solve(t(ans$rotmat) %*% ans$rotmat)
のようにしたときと，
factanal(x, factors=n, rotation="promax")
で，答えが違う。

「順列検定」って何? その2

「R による計算機統計学」が発売されていた。

「順列検定」は「並べ替え検定」に，

「確率変数を生み出す方法」は「確率変数の発生方法」に直っていた。

少し手を加えておかないとね

裏で何がどのようにやられていようと，一つの関数を呼ぶだけで描画したいとでも?以下のような関数を定義しておいて，plot(glmオブジェクト) とやれば，グラフは描かれます。 -- 河童の屁は，河童にあらず，屁である。? 2011-06-03 (金) 14:18:45

plot.glm <- function(a)
{
    x <- a$data$x
    y <- a$data$y
    plot(x, y)
    points(x, a$fitted.values, col="red")
    d <- data.frame(x=seq(x[1], x[length(x)], length=500))
    lines(d$x, predict(a, d, type="response"), col="red")
}
# 使用例
d <- data.frame(x=1:10, y=c(0,0,0,1,0,1,1,1,1,1))
a <- glm(y~x, data=d, family=binomial)
plot(a)
# 当然ながら，logit 以外にも対応
plot(glm(y~x, data=d, family=binomial(link="probit")))
plot(glm(y~x, data=d, family=gaussian))

は，関数外部で定義されるデータフレームの変数名が x, y でないときには正しく動かないっぽいので，取りあえず以下のように変更しておかないと行けないのではないかなあ。
でも，eval(parse(text="...")) を使うのはちょっと（かなり）ださいなあ。

plot.glm <- function(a)
{
    c.x <- as.character(a$term[[3]])
    c.y <- as.character(a$term[[2]])
    x <- eval(parse(text=sprintf("a$data$%s", c.x)))
    y <- eval(parse(text=sprintf("a$data$%s", c.y)))
    eval(parse(text=sprintf("plot(x, y, xlab='%s', ylab='%s')", c.x, c.y)))
    points(x, a$fitted.values, col="red")
    d <- data.frame(x=seq(x[1], x[length(x)], length=500))
    eval(parse(text=sprintf("colnames(d) <- '%s'", c.x)))
    lines(d[,1], predict(a, d, type="response"), col="red")
}
d <- data.frame(a1=1:10, a2=c(0,0,0,1,0,1,1,1,1,1))
a <- glm(a2~a1, data=d, family=binomial)
plot(a)

横着するな

ロジスティック回帰の曲線のグラフを描きたい †

うー? (2011-06-02 (木) 17:03:52)

ロジスティック回帰をしています。
量xと反応の割合pについて、glm関数でロジスティック回帰式を出しました。
xとpについてプロットを行い、この図に回帰曲線を加えたいと思います。

単回帰であればabline関数とlm関数を用いてコマンド一つで回帰直線を描くことができますが、ロジスティック回帰についてはそれが見当たりません。

いくつか本やウェブページを調べたのですが、独立変数を非常に細かくとってpredictとpointsで作図するとか、glmで求めた予測値についてlines関数で補間するなど、ちょっと回りくどい感じがします。

curve関数で書いてしまおうかとも思ったのですが、高水準作図関数なので実データとの作図の順番が逆になってしまいます。

何かいい方法はないでしょうか？　自分で関数を書く以外ないでしょうか？

 

curve だって，内部では seq.int なんか使っているんだから，同じでしょう。作図の順番が逆って，そんなの関係ないでしょう。できあがったものが，どんな順序で書かれたかなんか無関係。

要するに，少しの面倒もしたくないというわがままを言っているのか。