2010年11月10日

作図時に任意の数字や文字に使用する軸を変更する方法

後輩からよくわからないという連絡を受けたので、自分でも再確認のためにメモ。
確かに検索しても見つけにくいような気もする。
plot(seq(1, 10), xaxt = "n", yaxt="n")
axis(1, at=1:10, LETTERS[1:10])
axis(2, at=1:10, seq(100, 1000, length=10))

出来た図
change_axis.pdf

ラベル:R
posted by しばきん at 22:16| 神奈川 ☀| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする

2010年11月02日

多重共線性の指標VIF

多重共線性の指標でVIFというものがあります。例えば青木先生のサイトに載っています。
http://aoki2.si.gunma-u.ac.jp/R/tolerance.html
相関係数を利用して計算します。10以上だと多重共線性の指標だそうです(ちなみにそれを証明している文献は見たことがありません)。相関係数が高いから一緒に回帰式に含めるべきではない、とか言う人がいるが、実際にはどのくらいの相関係数が高い、ということになるのだろう。ちょっとやってみました。

> set.seed(1)
> a <- rnorm(10)
> b <- rnorm(10)
> c <- rnorm(10)
> res <- lm(a ~ b + c)
>
> ###単相関を計算###
> co <- cor(b, c)
>
> ###VIFを計算###
> library(DAAG)
> vif(res)
b c
1.5744 1.5744
>
> ###確認###
> 1/(1-co^2)
[1] 1.574429
>
> ###要は10=1/(1-x^2)を解けば良いのだから解くとx^2=9/10###
> sqrt(9/10)
[1] 0.9486833


意外にも相関係数0.948までは共線性がないことになるらしい。
これは・・・高く感じるのだけどそうでもないのかなぁ?
ちなみに変数を2つ以上に増やすと当然ながら単相関の結果とは合わなくなります。

> set.seed(1)
> b <- seq(1, 10)
> c <- seq(1, 10) + rnorm(10, 0, 1)
> d <- seq(1, 10) + rnorm(10, 0, 1)
> res <- lm(a ~ b + c + d)

> vif(res)
b c d
45.927 22.478 12.110

> cor(b, c)
[1] 0.9726364
> cor(b, d)
[1] 0.9485793
> cor(c, d)
[1] 0.8917848

ま、ちゃんと逆行列から計算しろってことですかな。こんな感じ?

> cor_mat <- matrix(c(1, cor(b, c), cor(b, d),
+ cor(c, b), 1, cor(c, d),
+ cor(d, b), cor(d, c), 1),
+ 3, 3)

###自分で3*3逆行列を自力で解くのは嫌なので、勝手に計算させました###
> solve(cor_mat)[1, 1]
[1] 45.92692
> solve(cor_mat)[2, 2]
[1] 22.47831
> solve(cor_mat)[3, 3]
[1] 12.10953


vif(res)の結果と一緒になりました。青木先生のサイトの言葉を借りると、「独立変数 i を残りの独立変数で予測するときの重相関係数になっている」ということなので、単相関で0.948以下だから大丈夫、とは思わないようにしたいと思います。

#2011/11/16日追記
VIFというキーワードでこのサイトに来る方が多いようなので、補足です。
つまり相関係数だけ見ても、どれくらいの相関係数なら多重共線性が生じていると言えるのか、という解はないので、ちゃんとVIFを計算して判断しましょう、という結論です。
ラベル:R 統計学
posted by しばきん at 18:18| 神奈川 ☁| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2010年11月01日

サンプル数が多いと本当に有意差が出やすいのか

こんなタイトルをつけるあたり、検定の仕組みとか標準誤差を理解してないのでは、とか思われそうですが、まぁ何事も理論だけじゃなく自分でも確認ってことで。分散は1で固定して、平均値0.01と0.02の違いをちゃんと検出出来るのかやってみました。

set.seed(12)

iter <- 100

###サンプル数###
c <- c(10, 1000, 100000)

d <- matrix(NA, iter, length(c))

for (j in 1 : length(c)) {
 for (i in 1 : iter) {
  a <- rnorm(c[j], 0.01, 1)
  b <- rnorm(c[j], 0.02, 1)
  
###t検定してみる###
  d[i, j] <- t.test(a, b)$p.value
 }
}

apply(d, 2, mean)

#pdf("Num_P.pdf")
par(mfrow=c(3, 1))
for(i in 1 : length(c)) {
hist(d[, i], breaks=20, col=c(2, rep(0, 19)), xlab="p_value", main=
eval(
parse(
text=
paste(c("c[", i, "]", sep="")
)
)
)
)
}
#dev.off()

出来た図。赤になっているところが有意差がちゃんと出たとこです。サンプル数が多いと確かに検出できた割合があがるようです。
Num_P.pdf
ラベル:R 統計学
posted by しばきん at 11:46| 神奈川 ☔| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。