StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

Rのsource()関数で呼び出すスクリプトに引数を渡すとき

学生に説明する必要が発生したためエントリ起こしておく。
コマンドライン引数みたいな感じで、source()で呼び出されるスクリプトに呼ぶ側のスクリプトから引数を渡すときは、以下のようにやればよい。
source1.Rからsource3.Rまでのスクリプトを準備してあり、それをreader.Rからsource()関数で呼び出す。
その時に、commandArgs()という関数を使う。

#呼ばれる側のスクリプト1(source1.R)
# 文字列が1つ渡される想定

print(commandArgs())
#呼ばれる側のスクリプト2(source2.R)
# ベクトルが一つ渡される想定

mean(commandArgs())
#呼ばれる側のスクリプト3(source3.R)
# 文字列と数値がリストで渡される想定

args <- commandArgs()
print(rep(args[[1]], args[[2]]))
# 呼ぶ側のスクリプト(reader.R)

# 渡す引数が1個なら単にそれを書けばいい
commandArgs <- function(...) {'ばか'}
source('source1.R')

# ベクトルを渡すこともできる
commandArgs <- function(...) {c(1,2,5,7,8)}
source('source2.R')

# 複数の引数をリストで渡す
commandArgs <- function(...) {list('あほ', 5)}
source('source3.R')


なおここで、function(...)の三点ドットは、定義されていない引数を表す記号。
以下、実行結果。

> # 渡す引数が1個なら単にそれを書けばいい
> commandArgs <- function(...) {'ばか'}
> source('source1.R')
[1] "ばか"
> 
> # ベクトルを渡すこともできる
> commandArgs <- function(...) {c(1,2,5,7,8)}
> source('source2.R')
[1] 4.6
> 
> # 複数の引数をリストで渡す
> commandArgs <- function(...) {list('あほ', 5)}
> source('source3.R')
[1] "あほ" "あほ" "あほ" "あほ" "あほ"

researchmapのcsv

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)みたいに"引数名="をつけるとエラーになる。

相関係数の差の検定と、回帰係数の差の検定を、Rでやる

たまに、2つの相関係数が有意に異なるのかや、1つの重回帰モデル中の2つの回帰係数が有意に異なるかを示せると、主張が通りやすいという場面がある。
まぁ、あまり必要になることはないのだが、相関係数の差の検定や回帰係数の差の検定について、日本語でググると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ダミーになってるという意味)で、有意になってるので、さっきの回帰係数の差は有意ってことになります。

*1:何かで読んだけど出典を忘れてしまった。

*2:回帰係数の差は交互作用で検定すればいいよっていうやり方自体は、以前、SPSSのマニュアルかなにかで読んだのだが、出典は忘れてしまった。でもまあ心理統計では共分散分析とかの解説でよく出てくる内容だと思います。

「不偏分散の平方根」は不偏標準偏差ではない件(メモ)

不偏分散の平方根を取っても不偏標準偏差にはならないという話があり、私は不真面目な研究者なのでそもそもそんなこと考えたこともなかったですが、知り合いが「数学的な導出はみれば分かるが、“平方根を取っては駄目な理由”が直観的に理解できなくて気持ち悪い」と言っていました。


ちなみに↓の記事によると、「不偏分散の平方根」を「不偏標準偏差」と記述しているケースはけっこうあるらしく、統計ソフトの計算もそうなってるらしいです(自分では確認してないですが)。
不偏分散の平方根は不偏標準偏差じゃなかった - 静粛に、只今統計勉強中
不偏標準偏差とは?:統計検定を理解せずに使っている人のために


で、ググってみたら下記のような解説があり、
Econometrics Beat: Dave Giles' Blog: Unbiased Estimation of a Standard Deviation
要するに直観的には、

  • 不偏推定というのは「期待値」に一致する値を求める手続きである
  • 期待値には線形性がある(たとえばE(X+Y)=E(X) + E(Y)
  • 平方根を取るのは非線形な変換である
  • だから「不偏分散の平方根」は「母分散の平方根の期待値」にはならない

というような話みたいです。
15分ぐらいしか考えてないのでまちがってたらすいません。


それと合わせて思ったのは、不偏推定値というのはそもそも、母集団からランダムサンプリングして推定値を得るという手続きを何回も繰り返して平均を取ると、だんだん母集団の真の値に近づいていくという性質の推定値で、1回1回得られた「不偏標準偏差」や「不偏分散」というものは、「何かの標準偏差」「何かの分散」ではないと考えておくべきな気もする。95%信頼区間というのが、1回1回の推定で得られる「その区間」に何かの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))
}


使うと、こんな感じで昔のパソコンみたいな動きになる。


f:id:midnightseminar:20200920130355g:plain


ちなみに、環境によっては、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で指摘をもらって気づいたけど、↓の例だと縦方向にも処理できるのでベクトル化容易やなww 横方向にしか処理できない内容であとで試そう。

> ### 値を置換する処理を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