researchmapのcsv取り込みにクセがありすぎる。
アクションがinsertなのに"invalid_delete_reason,,削除理由が無効です。"というエラーが出るのはバグ？
ちなみに削除理由にmineとか設定すると通る。意味がわからん。
そもそも、csvなのにヘッダーの上に"published papers"とか書かせるフォーマットなのが狂気。
Stanでよく忘れる、よく間違える書き方
- 文末の;の忘れ。
- ファイルの最後は空の行に
- forループの範囲の1:NのところのNを、intじゃなくてrealで宣言してしまっている。
- transformed parametersのブロックでforループを複数かくとコンパイルエラーになる理由がわからないときある。1つのループにまとめると通る。
- 定義の順序を間違えててまだ存在しない変数を参照してる。
- 自作関数を作って、つかうときに、引数の指定をRみたいにmy_function(x=1)みたいに"引数名="をつけるとエラーになる。
- rstanで"Error in serialize"とかのエラーがでるときは、一回Rを落としてプロセスをキルしてからやるといけたりする。
相関係数の差の検定と、回帰係数の差の検定を、Rでやる
たまに、２つの相関係数が有意に異なるのかや、１つの重回帰モデル中の２つの回帰係数が有意に異なるかを示せると、主張が通りやすいという場面がある。
まぁ、あまり必要になることはないのだが、相関係数の差の検定や回帰係数の差の検定について、日本語でググるとRでパッと使える方法がまとまっているわけではなさそうだったので、ここに書いておこう。
相関係数の差の検定
相関係数の差の検定は、Rの{psych}パッケージのr.test()関数で簡単にできる。
下記に関数の解説があるが、
https://www.rdocumentation.org/packages/psych/versions/2.0.12/topics/r.testr.test function - RDocumentation
Test4のtable1が崩れていたり、「Case A: where Masculinity at time 1 (M1) correlates with Verbal Ability .5 (r12), femininity at time 1 (F1) correlates with Verbal ability r13 =.4」とあるところはM1とF1が逆になっちゃってたり、出力の説明が不足していたりするので要注意。そこから参照されている論文Steiger(1980)をみると、いくつかの疑問は解決する。
1) この関数は1つの相関係数が有意かどうか検定したい場合にも使え、たとえば変数1・2ともにサンプルサイズが100で相関係数が0.66だとすると、
library(psych) r.test(n=100, r12=0.66)
というふうにする。多くの場合両側で検定すると思うが、片側で検定したい場合は
twotailed = F
という引数を付ければよい。
2) ここから2つの相関係数の差の話に移る。まずは、2つの独立した（サンプルを共有していない）相関係数の差を検定する場合。サンプルサイズは異なっても良い。
変数1と2はサンプルサイズが100で相関係数が0.3、変数3と4はサンプルサイズが120で相関係数が0.25だった場合、
r.test(n=100, r12=0.30, r34=0.25, n2=120)
サンプルサイズが等しいのであれば、n2は付けなくてもよい（デフォルトでnと同じになる）。
3) サンプルを共有している変数が3つ（1・2・3）あって、1と2、1と3というふうに1つの変数を共有する形で相関係数同士の差を検定したい場合。
r.test(n, r12, r23, r13)
というふうにする。この場合、r12とr13を比べた検定結果が返される。このとき、r23も引数として与える点に注意。（使い方の説明が不親切だが、冒頭の解説リンクのtest4のCaseAと同じことで、そちらについては論文Steiger(1980)と照合すると分かる。）
4) サンプルは共有してるけど変数は共有していないような、2組の相関係数を比べる場合。（冒頭の解説リンクでいうとtest4のCaseBに相当）
r.test(n, r12, r34, r23, r13, r14, r24)
冒頭の解説リンクは使い方の説明が不親切なのだが、論文Steiger(1980)と照合すると分かる。差を検定したい関心のある相関係数がr12とr34で、この2つの相関係数の差の検定結果を返している。
他に参考になるページ
http://www.snap-tck.com/room04/c01/stat/stat05/stat0501.html統計学入門−第5章
母相関係数の差の検定の概要と結果 :: 【公式】株式会社アイスタット｜統計分析研究所
回帰係数の差の検定
重回帰分析の回帰係数の差の検定については、たとえば複数の回帰係数の信頼区間を出して、大きい方の下限値と小さい方の上限値が重なるかどうかという基準で検定すると、（たとえば5%水準の検定をしたいときに95%信頼区間を用いると）厳しすぎる検定になるらしい。*1
で、使えるケースが多少限られるものの、心理統計の分野で分散分析の延長で共分散分析を習うときに出てくる「平行性の検定」の考え方で、「交互作用が有意になるかどうか」という観点で検定すると簡単にできる。*2
d <- airquality # 練習用データセットの読み込み d <- d[complete.cases(d),] # 欠損値削除 summary(lm(Ozone~Wind+Temp, data = d))
このようにすると以下のような結果が得られる。
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -67.3220 23.6210 -2.850 0.00524 ** Wind -3.2948 0.6711 -4.909 0.0000032617283 *** Temp 1.8276 0.2506 7.294 0.0000000000529 ***
この-3.29と1.83が有意に異なるかを知りたいという話で、まあ標準誤差とかをみても余裕で異なりそうではあるが、以下のように（従属変数であるOzoneは残して）ロング型にデータを変更し、交互作用の検定をすればいい。これは要するに、2つの変数を1つの変数にまとめて、変数名を新たな変数（ダミー変数）として設け、いわば「値」と「変数名」の交互作用をみるような感じ。
d2 <- d %>% select(Ozone, Wind, Temp) %>% # 変数を限定 pivot_longer(-Ozone, names_to = 'Variable', values_to = 'Value') # ロング型に変換 summary(lm(Ozone~Variable*Value, data = d2))
以下のような結果が得られる。
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -147.6461 19.7613 -7.471 0.0000000000019 *** VariableWind 246.6873 21.0072 11.743 < 0.0000000000000002 *** Value 2.4391 0.2522 9.673 < 0.0000000000000002 *** VariableWind:Value -8.1679 0.7210 -11.329 < 0.0000000000000002 ***
「VariableWind:Value」のところが交互作用（Tempが基準=0になって要するにWindダミーになってるという意味）で、有意になってるので、さっきの回帰係数の差は有意ってことになります。
後で知ったのですが、回帰係数の差の検定については、{car}パッケージのlinearHypothesis関数を使って、簡単に調べることができる。これは、2つの説明変数の回帰係数が等しいという帰無仮説を検定するもの。
あまりよく考えてないんですが、従属変数を複数もつ多変量回帰（Multivariate Regression。carパッケージのManovaで検定したりする）の場合、「すべてのモデルにおいて両変数の回帰係数に差がない」という帰無仮説を検定してるような気がする。モデルの数が増えると、帰無仮説がほとんど棄却されるようになります。上のように交互作用項を入れた上でManovaする方法の場合、交互作用はそれほど有意になりやすくはないです。
# 必要なパッケージの読み込み library(car) # airqualityデータセットを使用 data("airquality") # NA値の除去 airquality_complete <- na.omit(airquality) cor(airquality_complete) # 線形回帰モデルの作成 model <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality_complete) # Solar.R と Wind の回帰係数が等しいという仮説を検定 linearHypothesis(model, "Solar.R = Wind") # Multivariate Regression（多変量回帰：従属変数が複数ある）でもできる model_multi <- lm(cbind(Ozone, Wind) ~ Solar.R + Temp + Month + Day, data = airquality_complete) # Manovaで説明変数の総合的な有意性を確認 Manova(model_multi) # Solar.R と Tempの回帰係数が等しいという仮説を検定（棄却できる） linearHypothesis(model_multi, "Solar.R = Temp") # Month と Dayの回帰係数が等しいという仮説を検定（棄却できない） linearHypothesis(model_multi, "Month = Day")
「不偏分散の平方根」は不偏標準偏差ではない件（メモ）
不偏分散の平方根を取っても不偏標準偏差にはならないという話があり、私は不真面目な研究者なのでそもそもそんなこと考えたこともなかったですが、知り合いが「数学的な導出はみれば分かるが、“平方根を取っては駄目な理由”が直観的に理解できなくて気持ち悪い」と言っていました。
ちなみに↓の記事によると、「不偏分散の平方根」を「不偏標準偏差」と記述しているケースはけっこうあるらしく、統計ソフトの計算もそうなってるらしいです（自分では確認してないですが）。
不偏分散の平方根は不偏標準偏差じゃなかった - 静粛に、只今統計勉強中
不偏標準偏差とは?：統計検定を理解せずに使っている人のために
で、ググってみたら下記のような解説があり、
Econometrics Beat: Dave Giles' Blog: Unbiased Estimation of a Standard Deviation
要するに直観的には、
- 不偏推定というのは「期待値」に一致する値を求める手続きである
- 期待値には線形性がある（たとえばE(X+Y)=E(X) + E(Y)）
- 平方根を取るのは非線形な変換である
- だから「不偏分散の平方根」は「母分散の平方根の期待値」にはならない
というような話みたいです。
15分ぐらいしか考えてないのでまちがってたらすいません。
それと合わせて思ったのは、不偏推定値というのはそもそも、母集団からランダムサンプリングして推定値を得るという手続きを何回も繰り返して平均を取ると、だんだん母集団の真の値に近づいていくという性質の推定値で、１回１回得られた「不偏標準偏差」や「不偏分散」というものは、「何かの標準偏差」「何かの分散」ではないと考えておくべきな気もする。95%信頼区間というのが、１回１回の推定で得られる「その区間」に何かの95%が収まるという意味での区間ではないのと同様。
Rのループ中に進捗率を表示するプログレスバーを作る
Rのループで使えるプログレスバーは、いくつかのパッケージで提供されているみたいなのですが、自分で書くのも簡単なので、単純な関数でつくってみた。
進捗が知りたいのは時間のかかる処理をするときであり、時間がかかるなら無駄な計算は省きたいので、ループ1回ごとには出さずに100回とか1000回ごとに表示したほうがいいかもしれない。
### 自作プログレスバーをつくる # バーは40文字からなり、-が#に置き換わっていくようにした。 # 右端に％表示を出しておく。 # 処理は単に、ループ全体の何番目をやってるのかをwhichとlengthで # 取得して、それを文字列でのバーに置き換えてmessageとして吐くだけ。 # '\r'は行頭復帰を意味し、appendLF=FALSEは改行しないことを意味する。 # これらの組み合わせで、「表示をまるごと更新する」の意味になる。 # '\r'は、message()内の引数としては、prg.barの前でも後ろでもいいし、 # paster時にくっつけといてもよい。 # （messageは、複数の文字列を任意個数並べられる。） show.progress <- function(i, x){ # iとxにはforのループ変数とリストを与える prg <- round(which(x==i)/length(x)*100) done <- paste(rep('#', round(prg/2.5)), collapse='') remain <- paste(rep('-', 40-round(prg/2.5)), collapse='') prg.bar <- paste('|', done, remain, '| ', as.character(prg), '%', sep='') message('\r', prg.bar, appendLF=FALSE) } ### 試用してみる # 素因数分解できるパッケージ library(gmp) factorize(as.bigz("5656")) # 1億〜1億3万までの素因数分解を進捗バーつきでやってみる x<- 100000000:100030000 primes <- list() for (xi in x){ show.progress(xi,x) xi <- as.character(xi) p <- factorize(as.bigz(xi)) p <- as.numeric(p) primes <- c(primes, list(p)) }
使うと、こんな感じで昔のパソコンみたいな動きになる。
ちなみに、環境によっては、messageの行のあとにflush.console()を挟まないと、100%まで行ってからいきなり表示されるような場合もあるかもしれない。自分の場合は不要だったが。
for、apply、ベクトル演算の処理速度の比較
ツイッターのライムラインで、forループをapplyに置き換えた場合の高速化の話が流れていて（こちら）、気になって検索したところ、applyよりむしろwithを使えと言っている人がいた。
r - apply() is slow - how to make it faster or what are my alternatives? - Stack Overflow
ちなみにこれは、with関数が早いというより、要するに「データフレームの行ごとの計算」を「ベクトル同士の計算」に置き換えることによって速くなっていると考えたほうがよい。withのメリットは、データフレームの列を指定するときに'd$'を書かなくでよくなるという点にある。
以下では、1万行×3列のデータで、横向きに標本分散を計算して1万個取り出す処理を、
- forループで空のベクトルにアペンドしていく
- forループで長さを指定したベクトルを更新していく（Pythonのリストの処理でこうすると速くなる場合があったのでやってみる）
- applyを使う
- applyを使い、かつ関数を事前にコンパイルしておく
- ベクトル同士の計算にする（$で列を取り出す）
- ベクトル同士の計算にする（indexで列を取り出す）
- ベクトル同士の計算にする（with関数をつかう）
の7通りで試してみた。
# パッケージよみこみ library(rbenchmark) # 速度を測る作業のラッパー library(compiler) # 関数をコンパイルする # データフレームさくせい d <- data.frame( x1 = runif(10000), x2 = runif(10000), x3 = runif(10000) ) # 標本分散を出す関数を定義しておく svar <- function(x){ return(sum((x-mean(x))^2)/(length(x))) } # 関数をコンパイルする c.svar <- cmpfun(svar) # 処理時間を比較する bm1 <- benchmark( 'for_append' = { v1 <- c() for(i in 1:nrow(d)){ v1 <- c(v1, svar(unlist(d[i,]))) } }, 'for_replace' = { v2 <- rep(NA, nrow(d)) for(i in 1:nrow(d)){ v2[i] <- svar(unlist(d[i,])) } }, 'apply' = { v3 <- apply(d, 1, svar) }, 'apply+compile' = { v4 <- apply(d, 1, c.svar) }, 'vectorize_$' = { v5 <- ((d$x1-(d$x1+d$x2+d$x3)/3)^2 + (d$x2-(d$x1+d$x2+d$x3)/3)^2 + (d$x3-(d$x1+d$x2+d$x3)/3)^2)/3 }, 'vectorize_index' = { v6 <- ((d[,1]-(d[,1]+d[,2]+d[,3])/3)^2 + (d[,2]-(d[,1]+d[,2]+d[,3])/3)^2 + (d[,3]-(d[,1]+d[,2]+d[,3])/3)^2)/3 }, 'vectorize_with' = { v7 <- with(d, ((x1-(x1+x2+x3)/3)^2 + (x2-(x1+x2+x3)/3)^2 + (x3-(x1+x2+x3)/3)^2)/3) }, replications = 100) # いらん列を消して経過時間でソートしておく bm1 <- bm1[1:4] %>% arrange(elapsed)
replications = 100というのは100回処理を繰り返して比較してるということ。
以下のような結果が得られた。
test replications elapsed relative 1 vectorize_with 100 0.032 1.000 2 vectorize_$ 100 0.034 1.063 3 vectorize_index 100 0.042 1.313 4 apply+compile 100 6.618 206.812 5 apply 100 7.565 236.406 6 for_replace 100 52.767 1648.969 7 for_append 100 93.641 2926.281
elapsedが経過時間を表し、relativeは最速のものとの比を示したもの。
とにかく、ベクトル同士の演算に置き換えられる場合はめちゃめちゃ速くなることが分かる。アクセスの仕方で差はほぼ無く、withを使うと記述が簡単になるのでwithを使っておけばよさそうだ。
applyはforループより1桁速いが、それでもベクトル演算に比べると2桁倍の時間がかかっている。コンパイルしても凄く速くはならない。
forループは、上書き方式にするとアペンド方式より速いが、いずれにしてもベクトル演算の数千倍の時間がかかっている。
↓のように値を置き換えるような処理だと、ベクトル演算にはできないので、とりあえずforよりはapplyを使っておくのが良さそうな気がする。コンパイルしたほうが遅くて草。【追記】Twitterで指摘をもらって気づいたけど、↓の例だと縦方向にも処理できるのでベクトル化容易やなｗｗ 横方向にしか処理できない内容であとで試そう。
> ### 値を置換する処理をforとapplyで比較する > > # ベクトル中の0.5より小さい値を0に置換する関数 > rpl <- function(x){ + replace(x, which(x < 0.5), 0) + } > > c.rpl <- cmpfun(rpl) > > bm2 <- benchmark( + 'for' = { + d.new1 <- d + for(i in 1:nrow(d)){ + d[i,] <- rpl(d[i,]) + } + }, + 'for+compile' = { + d.new2 <- d + for(i in 1:nrow(d)){ + d[i,] <- c.rpl(d[i,]) + } + }, + 'apply' = { + d.new3 <- d + apply(d.new2, 1, rpl) + }, + 'apply+compile' = { + d.new4 <- d + apply(d.new2, 1, c.rpl) + }, + replications = 5 + ) > > bm2[1:4] %>% arrange(elapsed) test replications elapsed relative 1 apply 5 0.246 1.000 2 apply+compile 5 0.247 1.004 3 for 5 20.957 85.191 4 for+compile 5 22.890 93.049
ggplot2で2軸グラフを描く時の軸スケーリングの作業
ggplot2で2軸のグラフを描くときは、先日のエントリでも書いたように、ggplot2自身は左軸（第1軸）と右軸（第2軸）を別々の情報として持つことはできないので、左軸と右軸の尺度の違いを自分で設定して変換しなければならない。
あとで使いまわすので、このスケーリングの作業を行うスケーラを以下のように定義しておいた。
最初は1行で書いてggplotの描画の中に埋め込んでたけどあとで混乱しないように分けて書いておいた。
library(ggplot2) data(airquality) # 練習用データのよみこみ ### 変数を追加 # x軸用に月と日の列を日付の変数にしておく（2行に分けてかいてる） airquality <- airquality %>% mutate(Date = paste(as.character(Month), '/', as.character(Day), sep='')) %>% mutate(Date = as.Date(Date, format='%m/%d')) ### y1, y2の目盛り範囲を決めておく # ここでは恣意的に決めてるが、最大と最小を取るとかでもいいと思う y1.lim <- c(0, 25) y2.lim <- c(50, 100) ### スケーラの関数を書いておく # 上でつくった、y1とy2の目盛りの範囲を定めた要素数2のベクトルを使って、 # いい感じにスケールを合わせる。 # 変数のスケーラ。 # pにy2の値ベクトルを与えると、y1の尺に合わせた数字に変換。 # y2をまずゼロ基準に戻し、y2とy1のlimの幅の比でスケーリング # した後で、y1のゼロ基準からの乖離分を足す。 variable_scaler <- function(p, lim1, lim2){ to_zero <- p-lim2[1] y1_range <- lim1[2]-lim1[1] y2_range <- lim2[2]-lim2[1] scaled <- to_zero*y1_range/y2_range from_zero <- scaled + lim1[1] return(from_zero) } # 第2軸の目盛りのスケーラ。 # pは、sec_axis()の'.'になる。y1の目盛りをy2の目盛りに読み替えるもの。 # y1の目盛りをまずゼロ基準に戻し、y1とy2のlimの幅の比スケーリング # した後で、y2の目盛りのゼロ基準からの乖離分を足す。 axis_scaler <- function(p, lim1, lim2){ to_zero <- p-lim1[1] y1_range <- lim1[2]-lim1[1] y2_range <- lim2[2]-lim2[1] scaled <- to_zero*y2_range/y1_range from_zero <- scaled + lim2[1] return(from_zero) } ### 描画してみる airquality %>% ggplot(aes(x=Date)) + geom_line(aes(y=Wind, colour='Wind')) + geom_line(aes(y=variable_scaler(Temp, y1.lim, y2.lim), colour='Temp')) + scale_y_continuous(limit=y1.lim, # 第1軸の範囲 breaks=c(0, 5, 10, 15, 20, 25), # 第1軸の目盛り sec.axis=sec_axis( ~(axis_scaler(., y1.lim, y2.lim)), # 軸スケーリング breaks=c(50, 60, 70, 80, 90,100), # 第2軸の目盛り name='Temp') # y2のラベルはここで設定する ) + labs(title = 'DUAL AXIS CHART', x='Date', y='Wind', colour = 'Variable')
もちろんy2のラベルもscale_y_continuousの中でnameをつかって指定してもよい。
2つ目の折れ線のyを指定するときに変数スケーラを使い、sec.axisを指定するときにformulaの中で軸スケーラを使う。
