単なるメモですが、
rm(list=subset(ls(), !(ls() %in% c('hoge', 'hage'))))
rm(list=ls()[-which(ls() %in% c('hoge', 'hage'))])
'hoge', 'hage'のところに、消したくないオブジェクト名を書いておく。
後者のほうが2文字短い。
Rでファイルをダウンロードするかどうか確認させる関数
ネット上にあるcsvファイルをダウンロードしてきて使う場合、read.table()にurlを与えて直接データフレームをつくってしまう場合もあれば、ファイルとしてダウンロードして置いておきたい場合もあります（実行のたびにダウンロードしたくない等の理由で）。
で、実行するときに、すでにダウンロード済みなのであれば再ダウンロード（上書き）はしたくないという場合があるので、確認しながらダウンロードできると便利です。
簡単な処理ではありますが、コードの使い回しのためにここに書いておきます。
そのファイルがすでに存在するかどうかをメッセージで表示した上で、ダウンロードするかどうか訊いてくるので、コンソールでyを入力したらダウンロード、nを入力したら中止ということにしておきます。
ファイルがすでに存在する場合、そのファイルの更新日時を確認してから（たとえば古かったら）ダウンロードしたいという場合もあるので、更新日時も取得して表示するようにしておきます。
処理確認ダイアログを出す関数askYesNo()の使い方は、askYesNoのドキュメントをみてください。
confirm_file_dl <- function(url,dst){ # dstには保存するときのパス＆ファイル名を与える if(file.exists(dst)) { message(paste('\n\"', dst, '\"', ' already exists!', sep='')) mtime <- file.info(dst)$mtime # ファイルの更新日時の取得 message(paste('(Last modified at ', mtime, ')\n', sep='')) } else { message(paste('\n\"', dst, '\"', ' doesn\'t exist!\n', sep='')) } dl_or_not <- askYesNo(msg='Do you want to download the file?', default = FALSE, prompts = 'y/n/na') if(dl_or_not==TRUE) { download.file(url, destfile = dst) } }
たとえばGoogleのモビリティ・レポートのファイルは、国別データは1ファイルにまとまってるのですが、現時点で600MBもあって、外出中にテザリングで作業してる時なんかは、必要ないならダウンロードしたくないですね。
ファイルがすでに存在する場合。
> confirm_file_dl(url='https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv', + dst='source/Global_Mobility_Report.csv') "source/Global_Mobility_Report.csv" already exists! (Last modified at 2021-08-17 18:13:17) Do you want to download the file? (y/n/na)
ファイルが存在しない場合。（nを入力してダウンロードを実行）
> confirm_file_dl(url='https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv', + dst='source/Global_Mobility_Report.csv') "source/Global_Mobility_Report.csv" doesn't exist! Do you want to download the file? (y/n/na) y trying URL 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv' Content type 'text/csv' length 637979115 bytes (608.4 MB) ================================================== downloaded 608.4 MB
Rのsource()関数で呼び出すスクリプトに引数を渡すとき
学生に説明する必要が発生したためエントリ起こしておく。
コマンドライン引数みたいな感じで、source()で呼び出されるスクリプトに呼ぶ側のスクリプトから引数を渡すときは、以下のようにやればよい。
source1.Rからsource3.Rまでのスクリプトを準備してあり、それをreader.Rからsource()関数で呼び出す。
その時に、commandArgs()という関数を使う。
#呼ばれる側のスクリプト１（source1.R） # 文字列が1つ渡される想定 print(commandArgs())
#呼ばれる側のスクリプト２（source2.R） # ベクトルが一つ渡される想定 mean(commandArgs())
#呼ばれる側のスクリプト３（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)みたいに"引数名="をつけるとエラーになる。
- 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%が収まるという意味での区間ではないのと同様。
