R markdownを使ってみる

＃追記：どういうわけか、一度アップしてから、改行なり文字なり１文字でも追記すると画像の表示がうまくいくようだ。理由はわからない。

R markdownは知っていたが、使う必要性を感じなかった。しかし、blogをはじめた今となっては、使ってみても良いかも。

とりあえずやってみる

#以下markdownでの出力

マン・ホイットニーのU検定＝ウィルコクソンの順位和検定は、中央値の違いを見てるんじゃないよってことを研修医に軽く教えた時に使ったもの。

分散だけが違い、中央値や平均値には差がない場合をシミュレートする。

# 平均は同じで分散が違う
n <- 100  #n　cases
m <- 10000  #m times

d <- numeric(n)
dd <- numeric(n)
ttt <- numeric(n)
pp <- numeric(n)
wp <- numeric(n)
ptw <- numeric(n)
y <- x <- matrix(0, n * m, nrow = m, ncol = n)

x <- matrix(rnorm(n * m, 20, 1), nrow = m, ncol = n)
y <- matrix(rnorm(n * m, 20, 20), nrow = m, ncol = n)

for (i in 1:m) {
    d[i] <- mean(x[i, ]) - mean(y[i, ])
    dd[i] <- median(x[i, ]) - median(y[i, ])
    wt <- wilcox.test(x[i, ], y[i, ])
    wp[i] <- wt$p.value
    tt <- t.test(x[i, ], y[i, ], var.equal = T)
    ttw <- t.test(x[i, ], y[i, ], var.equal = F)
    pp[i] <- tt$p.value
    ptw[i] <- ttw$p.value
}

＃ｔ検定のP値の分布

hist(pp)

plot of chunk unnamed-chunk-2

length(pp[pp < 0.05])/m  #t test p

## [1] 0.0503

plot(sort(pp), xlab = "N", ylab = "p value of t test", ylim = c(0, 1))
abline(0, 1/m, col = 2)

plot of chunk unnamed-chunk-2

分散は違うがt検定はまずまず

＃マン・ホイットニーのU検定＝ウィルコクソンの順位和検定の分布

hist(wp)

plot of chunk unnamed-chunk-3

plot(sort(wp), xlab = "N", ylab = "p value of U test", ylim = c(0, 1))
length(wp[wp < 0.05])/m  #wilcox.test

## [1] 0.1008

abline(0, 1/m, col = 2)

plot of chunk unnamed-chunk-3

これはそれなりに有意差が出る。中央値に差はないがP<0.05になる例が１割位ある

つまり、マン・ホイットニーのU検定＝ウィルコクソンの順位和検定では、中央値に差がなくても有意差が出る。

ウェルチのt検定のP値の分布

hist(ptw)

plot of chunk unnamed-chunk-4

plot(sort(ptw), xlab = "N", ylab = "p value of welch t test", ylim = c(0, 1))
length(ptw[ptw < 0.05])/m  #t welch P

## [1] 0.0485

abline(0, 1/m, col = 2)

plot of chunk unnamed-chunk-4

これは良好。

まとめると、平均は同じ（中央値も同じ）だが、分散だけが違う二群だと

t検定やt検定welchの方法はP値が一様分布であり、P<0.05に

なるのは、５％であった。

しかし、マン・ホイットニーのU検定＝ウィルコクソンの順位和検定

だと、P値は明らかに小さい方に偏り、１０％近くの例でP<0.05となった。

なんか後ろにゴミが出る。良くわからないorz

<body>の間だけコピペしたら解決した

なんとか出力はできた。意外ときれいな記事になるようだ。

PCからは図が見えるが、スマホからだと見えない。

何かがおかしい。時間ができたら調査しよう。

r-statistics-fanの日記

統計好き人間の覚書のようなもの

R markdownを使ってみた（やり方は次回記事）

R markdownを使ってみる