2020年2月29日 (土)
Rをつかっているとformulaを書かないといけないことがあるけど(回帰系の関数を呼ぶときとかに)、なんか便利な書き方がありそうなのに全然使っていない。このたび、ふとstats::formula()のドキュメントを眺めていたら、恥ずかしながら「へぇー」と思うことがいくつかあった。
この話題、Rを使っている人にとってはたぶん馴染み深い話で、webをちょっと検索するだけで、丁寧に解説しておられるブログ記事をいくつもみつけることができる(たとえばこちら：[R] 予測モデルを作るには formula を活用せよ)。そういうのをきちんと読んできちんと理解しておけばいいんですけどね。そのときは「へぇー」と思うんだけど、すぐに忘れちゃうのである。
というわけで、今日は忘れないようにメモしておこう。なぜか執事風の丁寧語で。お嬢様、羊肉のソテーでございます、って感じで。
- 「
y ~ model」: 反応yをmodelで指定した予測子でモデル化いたします。 - 「
a + b」:aとbの線形和でございます。 - 「
a : b」:aとbの交互作用でございます。 - 「
a * b」:aとbのクロス、すなわちa+b+a:bでございます。 - 「
(a + b + c)^2」:(a + b + c)の二次のクロス、すなわちa,b,cの主効果と二次交互作用でございます [これ知らなかった... 文字通り展開するとa*a, b*b, c*cが生まれるけど、それらは算術的な意味での$a^2, b^2, c^2$ではなくて、単にa, b, cの主効果ってことになるのね。へー] - 「
a + b %in% a」:bはaにネストしております。すなわちa + a:bでございます。[ここちょっと不思議なところで...%in%は常に:と等価なのかしらん？] 「(a + b + c)^2 - a:b」:(a + b + c)^2からa:bを取り除いたもの、すなわちa + b + c + b:c + a:cでございます。- 「
y ~ x - 1」: 切片項を取り除いております。y ~ x + 0,y ~ 0 + xでもよろしゅうございます [←Rを使い始めたときから、formulaの右辺で0と1が等価だってのが気持ち悪くてしかたない。責任者でてこい] 「log(y) ~ a + log(x)」: このように、算術表現もお使いになれます [←私、formulaのなかでこういう変数変換を掛けるのが生理的に許せなくて、必ず自力で前処理している。なぜ気持ち悪いのかいま気が付いた。formulaのなかの演算子は算術演算子ではないのに(たとえば+は和ではないのに)、関数は算術表現として扱われるというのが気持ち悪いんだ...]- 「
a + I(b+c)」: aとb+cの線形和でございます。I()のなかの演算子は算術演算子として扱われるのでございます。 - 「
.」: 「そのformulaのなかにないすべての列」を表します [←formulaをupdateするときはちがう意味になるのだが、まあupdateなんて使わないからいいや] - 「
a + b + offset(z)」: aとbの線形和とzとの和でございます。zの係数は1となります。なお、この表現は受け付けない関数もございます。
ところで、lm(y ~ a/b)なんていう書き方をみたことがあるけど(y ~ a + a:bと等価だと思う)、stats::formula()のドキュメントにも、stats::lm()のドキュメントにも載ってなかった。あれえ？
雑記：データ解析 - 覚え書き：Rのformulaのなかで使える表現
仕事の都合でときどき混合効果モデルを組みたくなることがある。本腰を入れる場合はMplusを使うけど、ほんのちょっとした用事の場合など、ちょっぴり億劫に感じることもある。RStudioのなかでちゃちゃっと済ませちゃいたいな、なあんて... ううむ、Mplus信者としては敬虔さが足りない、すみませんすみません。
Rでなにかのモデルを組むとき、「どのパッケージを使えばいいんだ...」と困ることがあるけど、混合効果モデルはその典型例だと思う。古株としてはnlme (たしか標準パッケージだ)、lme4があり、さらにはglmmMLだの、glmmだの(あ、2018年で開発が止まったようだ)、かの複雑怪奇なMCMCglmmとか... 参っちゃいますね。
Bates, D., Machler, M., Bolker, B.M., Walker, S.C. Fitting Linear Mixed-Effects Models Using lme4.
というわけで、緊急事態に対応するため、どれでもいいからちょっと土地勘をつけておきたいと思い、lme4パッケージのビニエットに目を通した。発行年が書いてないんだけど、これの古い版は2015年にJ. Statistical Softwareに載っている模様。
ほんとはnlmeパッケージについて調べたかったのだけど、適切な文献が見当たらなかったのである。開発者による解説書、ピネイロ&ベイツ「S-PLUSによる混合効果モデル解析」を読めってことでしょうか。うーん、それは勘弁してほしい。あれは滅茶苦茶むずかしい。
ところが、今回はじめて知ったんだけど、lme4パッケージの筆頭開発者はDouglas Batesさん。ピネイロ&ベイツ本の第二著者、かつてピネイロさんとともにnlmeパッケージを開発した人であった。どうやら、Batesさんは2007年にnlme開発チームを離れ、ライバルとなる新パッケージlme4を作り始めた、ということらしい。へええ。
いわく。
lme4はnlmeと比べ...[さまざまな美点が挙げてある]。いっぽうnlmeは、残差に自己相関のような構造があるモデルとか、ランダム効果の共分散行列に構造があるモデルとかを組むときのUIが優れている。
まずは線形混合モデルから。
線形モデルは、ベクトル値ランダム反応変数$\mathcal{Y}$の分布
$\mathcal{Y} \sim N(\mathbf{X} \beta + \mathbf{o}, \sigma^2 \mathbf{W}^{-1})$
として記述できる。$\mathbf{W}$は対角行列で既知のウェイト、$\mathbf{o}$は既知のオフセット項。パラメータは$\beta$と$\sigma^2$である。
同様に、線形混合モデルは$\mathcal{Y}$とランダム効果ベクトル$\mathcal{B}$の分布として記述できる。$\mathcal{B}=\mathbf{b}$の下での条件つき分布は
$(\mathcal{Y} | \mathcal{B} = \mathbf{b}) = N(\mathbf{X} \beta + \mathbf{Zb} + \mathbf{o}, \sigma^2 \mathbf{W}^{-1})$
$\mathcal{B}$の分布は
$\mathcal{B} \sim N(\mathbf{0}, \mathbf{\Sigma})$
ただし$\mathbf{\Sigma}$は半正定値[すべての固有値が非負ってことね]。これは、分散成分パラメータ$\theta$に依存する相対共分散ファクター$\Lambda_\theta$を考えて
$\Sigma_\theta = \sigma^2 \Lambda_\theta \Lambda_\theta^\top$
と書くと便利である。[こういわれると抽象的で困っちゃいますが... ]
さて、lme4パッケージがご提供する線形混合モデルの関数は lmer() であります。
たとえば...[対象者の睡眠を9日間制限し、毎日反応時間課題をやったというデータの紹介...]。この場合、
lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
と書く。
実は、lmer() は以下の4つのモジュールの組み合わせである。
その1, Formula.
parsedFormula <- lFormula(formula = Reaction ~ Days + (Days | Subject), data = sleepstudy)
その2, 目的関数.
devianceFunction <- do.call(mkLmerDevfun, parsedFormula)
その3, 最適化.
optimizerOutput <- optimizeLmer(devianceFunction)
その4, 出力. (めんどくさいのでコードは省略するが、mkMerMod()という関数)
モジュール1, Formulaについて。
formula引数は
resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...
という形式をとる。FEexprは固定効果モデル行列$X$を決める表現。問題はその右側、ランダム効果の項である。
個々のREexprは線形モデルのformulaとして評価される。つまりlm()なんかのformulaと同じである。1と書いたら切片である。factorはRの因子型として評価される。縦棒は、効果とグルーピング因子との交互作用を表すと思えばよい。
いくつか例を挙げると、
(1 | g)は「gの個々の水準がランダム切片を持つ」という意味になる。ランダム切片の平均とSDがパラメータになる。1 + (1 | g)と書いてもよい。- ランダム切片の平均が既知の変数oであるならば
0 + offset(o) + (1|g)と書く。先頭の0は-1と書いてもよい。 - g1の下にg2があって、切片がg1の水準間でもg2の水準間でも動くなら
(1|g1) + (1|g1:g2)と書く。これは(1 | g1/g2)とも書ける。[←へー] - g1とg2がネストじゃなくてクロスしているなら
1 + (1|g1) + (1|g2)と書く。IRTなんかがそうだ(被験者と項目がクロスする)。先頭の1は略記できる。 - 切片と傾きの両方が動くなら
1 + x + (1+x|g)と書く。これはx + (x | g)と略記できる[←えっそうなの？ 知らなかった]。lme4では同じランダム効果項に出てくる係数はすべて相関することに注意。つまり切片と傾きは相関する。 - 切片と傾きを無相関にしないなら、
1 + x + (1|g) + (0+x|g)と書けばよい。これはx + (x || g)と略記できる[←へー！]。なお、こういう制約は予測子が比率尺度の場合に限るのが望ましい。この制約の下では予測子をシフトするだけで尤度が変わるから。[←ここ、ちょっとピンとこなかった... いずれ必要になったらゆっくり考えよう]
FEexprから$\mathbf{X}$, その右から$\mathbf{Z}$と$\Lambda_\theta$がつくられる。[以下、$\mathbf{Z}$と$\Lambda_\theta$の作り方について詳しい説明があったけど、半分くらいで力尽きた。メモ省略]
モジュール2, 目的関数モジュールは...[罰則付き最小二乗法によるあてはめの説明とかが12頁にわたって説明されている。読んだところで理解できそうにないのでまるごとパス]
モジュール3. 非線形最適化モジュールは... [ここもどうせわからんので読んでないけど、混合モデルに限らない汎用的な最適化関数らしい。ふーん]
モジュール4. 出力モジュールは... [単に出力の話かと思ったら、意外に仕事が多くて驚いた。まあこれは困ったときに読めばいいかな]
... 途中から面倒になってぱらぱらめくっただけなんだけど、まあいいや。
lme4パッケージの全貌についてのユーザ向け紹介を期待していたんだけど、このビニエットはlmer()の仕組みについて詳しーく説明するものであった。しょうがないのでメモしておくけど、lme4の主力関数は次の3つである。ときに、末尾のerってなんの略なんすかね？
- lmer(): 線形混合モデル
- glmer(): 一般化線形混合モデル
- nlmer(): 非線形混合モデル
nlmeパッケージでいうところのcorrelation = corAR1()のように、誤差項に自己相関を持たせる方法が知りたかったんだけどな... おそらく、lme4でもできなかないけど、めんどくさいのであろう。
論文：データ解析(2018-) - 読了：Bates, et al. (2015?) lme4::lmer()でたのしい線形混合モデリング
2020年2月17日 (月)
中国 人口減少の真実 (日経プレミアシリーズ)
[a]
村山 宏 / 日本経済新聞出版社 / 2020-01-11
皮肉でもなんでもなく素直に感心しているんだけど、日経新聞の編集委員ともなると、中国メディアの最近の面白記事をピックアップしてコメントをつけただけで、こうして本が一冊書けちゃうんですね。(本当に皮肉ではない。だって誰にでもできることではないじゃないですか)
ノンフィクション(2018-) - 読了：「中国 人口減少の真実」
わたしの名は赤〔新訳版〕 上 (ハヤカワepi文庫)
[a]
オルハン パムク / 早川書房 / 2012-01-25
わたしの名は赤〔新訳版〕 下 (ハヤカワepi文庫)
[a]
オルハン パムク / 早川書房 / 2012-01-25
年明けから身内の事情で、具体的になにができるわけでもないのにどうも気が休まらない日々が続いていたのだが、新幹線での移動中にずっと手に取っていたのが、以前から気になっていたオルハン・パムクの長編小説であった。良い本は少なからず心の支えとなってくれる。
いやあ、それにしても、かのノーベル文学賞作家の代表作がこんなエンタメ小説だとは。おそるべしパムク。
ポンス・ピラト ほか―カイヨワ幻想物語集
[a]
ロジェ カイヨワ / 景文館書店 / 2013-05-29
私鉄沿線の小さな書店でふと見かけて、R.カイヨワって「遊びと人間」の？ 社会学者じゃなかったの？ と意外の念に撃たれ、つい手に取った。
表題作は総督ピラトがイエスを処刑せよと迫られ、友人に相談したりするんだけど(友人はその後の西洋史を丸ごと幻視し、フランスの作家がこのやりとりを小説にするだろうとまでいう)、結局は処刑を拒否、おかげでキリスト教は興りませんでした...という人を食った小説であった。
ところで、この古風な体裁のペーパーバックを出した景文館書店なる版元は？ と疑問に思って検索してみたら、愛知県岡崎市にあるひとり出版社。立ち上げ早々に吉田知子の本を出版しようとして、作家に「あなたは財閥の御曹司かなにかか」と訝られたそうだ。ははは。
水の出会う場所
[a]
魚住陽子 / 駒草出版 / 2014-09-26
フィクション - 読了：「ポンス・ピラトほか」「水の出会う場所」「わたしの名は赤」
刑務所しか居場所がない人たち : 学校では教えてくれない、障害と犯罪の話
[a]
山本 譲司 / 大月書店 / 2018-05-17
「でも、断言しよう。[知的]障害のある人たちだって、かならずコミュニケーションができる。[...] いちばん共感しやすいのは「悲しい」という感情だ。障害のある人は、悲しい思いをたくさん経験しているから、他人の悲しみに敏感だ。[...] とくに、お母さんが悲しんでいる気持ちはすぐに察する。お母さんの表情が少しでも曇れば、彼らも落ち着かなくなる。そんなとき「悲しいね」「つらいね」って声をかけていると、いつしか心を開いてくれる」
コラボ゠対独協力者の粛清 (文庫クセジュ)
[a]
マルク・ベルジェール / 白水社 / 2019-11-09
「歴史認識」とは何か - 対立の構図を超えて (中公新書 2332)
[a]
大沼 保昭 / 中央公論新社 / 2015-07-24
話し手の先生(故人)は、村山内閣のときにできたアジア女性基金の理事を務めた人。私もまたあの頃は、国による補償ではなく募金だなんて、と反発していたのだが、このような世の中になってしまったいまにして思えば、なにもしないよりは全然ましだったんだろうなあ、と...
マトリ 厚労省麻薬取締官 (新潮新書)
[a]
瀬戸 晴海 / 新潮社 / 2020-01-16
マタイ受難曲 (ちくま学芸文庫)
[a]
礒山 雅 / 筑摩書房 / 2019-12-10
夜中にマタイ受難曲を流しながらちびちびと読んだ。楽曲解説にはいると、音楽学の専門用語が頻出しちんぷんかんぷんなんだけど、それでもそれなりに楽しく読了。
ノンフィクション(2018-) - 読了：「マタイ受難曲」「刑務所しか居場所のない人たち」「コラボ=対独協力者の粛正」「『歴史認識』とはなにか」「マトリ」
2020年2月16日 (日)
六四と一九八九:習近平帝国とどう向き合うのか
[a]
/ 白水社 / 2019-12-25
昨年に明治大で開かれたシンポジウムを基にした論文集だそうだ。張博樹という民主派知識人(現在は米在住) の章が興味深かった。
ハンナ・アーレント (ちくま新書)
[a]
大輔, 森分 / 筑摩書房 / 2019-06-06
ジャパン・ストーリー 昭和・平成の日本政治見聞録
[a]
ジェラルド・L・カーティス / 日経BP / 2019-05-23
新編 戦後翻訳風雲録 (大人の本棚)
[a]
宮田 昇 / みすず書房 / 2007-06-02
著者は早川書房の編集者出身。田村隆一がいかに迷惑な人であったかという話が縷々述べられていて、こういっちゃ申し訳ないけど、とても楽しく読んだ。
鮎川信夫が亡くなったあと、独身だったはずの鮎川の妻として年配の女性(最所フミ)が名乗りでて、周囲は仰天したのだそうだ。それって、ちくま文庫「英語類義語活用辞典」の最所フミではありませんか...
貞観政要 (ちくま学芸文庫)
[a]
呉 兢 / 筑摩書房 / 2015-09-09
ノンフィクション(2018-) - 読了：「六四と一九八九 習近平帝国とどう向き合うのか」「ハンナ・アーレント」「ジャパン・ストーリー」「新編 戦後翻訳風雲録」「貞観政要」
戦争は女の顔をしていない 1 (単行本コミックス)
[a]
小梅 けいと / KADOKAWA / 2020-01-27
Webで連載開始を知り仰天したんだけど、アレクシェービッチ「戦争は女の顔をしていない」のきわめて真摯なコミカライズである。
こうして単行本でまとめて読むと、圧倒されますね... きっと大変な評判を呼ぶだろう(もうとっくに呼んでいるのかもしれない)。
あさひなぐ (32) (ビッグコミックス)
[a]
こざき 亜衣 / 小学館 / 2020-01-30
きちじつごよみ 1 (フィールコミックス)
[a]
岩岡ヒサエ / 祥伝社 / 2020-01-24
あせとせっけん（６） (モーニングコミックス)
[a]
山田金鉄 / 講談社 / 2020-01-23
ハコヅメ~交番女子の逆襲~(11) (モーニング KC)
[a]
泰 三子 / 講談社 / 2020-01-23
いま一番気になるマンガ。これを読んでいると、マンガの巧さって一体なんなのか？としばし考え込んでしまう。著者は急激に腕を上げているとはいえ、いまだに「このコマはなんだかパースが狂っているのではないか」「人体というのはこういう風にはできていないのでは」などと思うことが多いんだけど、それでも引き込まれる。
コミックス(2020-) - 読了：「戦争は女の顔をしていない」「ハコヅメ」「あせとせっけん」「きちじつごよみ」「あさひなぐ」
アルテ ⑫ (ゼノンコミックス)
[a]
大久保圭 / 徳間書店 / 2020-01-20
北北西に曇と往け 4 (ハルタコミックス)
[a]
入江 亜季 / KADOKAWA / 2020-01-14
ゆりあ先生の赤い糸（５） (BE・LOVEコミックス)
[a]
入江喜和 / 講談社 / 2020-01-10
日々我人間２ (文春e-book)
[a]
桜 玉吉 / 文藝春秋 / 2020-01-16
東独にいた(1) (ヤンマガKCスペシャル)
[a]
宮下 暁 / 講談社 / 2019-12-20
ベルリンの壁崩壊の前夜を描く人間ドラマかと思いきや、超人的殺人能力者が暴れ回る少年誌的SFアクションであった...
コミックス(2020-) - 読了：「日々我人間」「東独にいた」「ゆりあ先生の赤い糸」「北北西に曇と往け」「アルテ」
めしばな刑事タチバナ 36 ガムの味 (トクマコミックス)
[a]
坂戸佐兵衛,旅井とり / 徳間書店 / 2019-12-27
おかあさんの扉9 謎の九歳児 (オレンジページムック)
[a]
伊藤 理佐 / オレンジページ / 2020-02-03
アンダーニンジャ（３） (ヤングマガジンコミックス)
[a]
花沢健吾 / 講談社 / 2020-02-06
重版出来! (14) (ビッグコミックス)
[a]
松田 奈緒子 / 小学館 / 2020-02-12
東京タラレバ娘 シーズン2(1) (KC KISS)
[a]
東村 アキコ / 講談社 / 2019-10-11
この著者のことだから、成算があってはじめた続編なんだろうけど...ううむ...まだちょっとよくわからない。
コミックス(2020-) - 読了：「重番出来！」「アンダー・ニンジャ」「東京タラレバ娘 シーズン2」「おかあさんの扉」「めしばな刑事タチバナ」
やはりこまめに記録した方がよさそうなので... 年明けから読んだ本、まずはコミックスから：
マイ・ブロークン・マリコ (BRIDGE COMICS)
[a]
平庫 ワカ / KADOKAWA / 2020-01-08
親友の自死を知った若い女が、友を苦しめた父親から遺骨を奪い、あてもなく鎮魂の旅に出る。いやあ、これはほんとに面白かった。亡き友が成人後は必ずしも良い友人ではなかった、という回想に胸を抉られる。
著者はこれが初単行本である由。絵柄には松本次郎さんの影響が非常に大きいと思うんだけど、それでもオリジナリティってのは顕れるものなんだなあ。
ペリリュー ―楽園のゲルニカ― 8 (ヤングアニマルコミックス)
[a]
武田 一義 / 白泉社 / 2020-01-29
パラオ・ペリリュー島での日本兵たちを描いた作品、舞台は終戦後の潜伏生活に移っている。
近所の最果て 澤江ポンプ短編集 (torch comics)
[a]
澤江 ポンプ / リイド社 / 2020-01-31
7人のシェイクスピア NON SANZ DROICT(11) (ヤンマガKCスペシャル)
[a]
ハロルド 作石 / 講談社 / 2020-02-06
さめない街の喫茶店
[a]
はしゃ / イースト・プレス / 2017-12-13
さめない街の喫茶店 2
[a]
はしゃ / イースト・プレス / 2019-11-17
コミックス(2020-) - 読了：「ペリリュー」「近所の最果て」「マイ・ブロークン・マリコ」「さめない町の喫茶店」「７人のシェイクスピア」
2020年2月14日 (金)
時系列データを扱っていると、たまーに単位根検定をやりたくなることがあるんだけど、Rのパッケージがいろいろあって困る。
仕方がないので、目についたやつについてメモを取った。なんというか、恥をさらしているような気がしますが...
tseriesパッケージ
時系列分析パッケージの古手だと思う。ADF検定, KPSS検定, PP検定の関数を持っている。
adf.test(x, alternative, k): H0「xは単位根を持つ」のADF検定。モデルは定数と線形トレンドを含む。引数は:
- x: 時系列
- alternative: 対立仮説の指定。値は"stationary", "explosive"。後者はnonstatonaryってこと？
- k: ラグ次数。デフォルトではtrunc((length(x)-1)^(1/3))となるそうだ。k=0とするとDickey-Fuller検定になる。
えーっと、いっつも混乱するんだけど、ここで「モデルは定数と線形トレンドを含む」っていうのはどういうことなの？
話を単純にするためにラグ次数を1として、ここで「定数と線形トレンドを含む」といっているのは、差分時系列について
$\Delta y_t = \beta_1 + \beta_2 t + (\phi-1) y_{t-1} + e_t$
というドリフト＋トレンドつきモデルを考えて$H_0: \phi = 1$の検定をいたします、ってこと？ それとも、元の時系列に切片と1次の項をいれて
$y_t = \mu + \beta_1 t + \phi y_{t-1} + e_t$
$\Delta y_t = \beta_1 + (\phi-1) y_{t-1} + e_t$
つまり差分時系列はドリフトつきモデルです、ってことなの？ うぐぐぐぐ。
[02/16 追記: コードをざっと眺めたところ、どうやら前者、つまり$\Delta y_t$のモデルに定数項と1次の項がはいるという話らしい... しらんけど]
kpss.test(x, null, lshort): H0「xはレベル定常ないしトレンド定常である」のKPSS検定。引数は
- x: 単変量時系列
- null: H0はレベル定常かトレンド定常か。値は"Level", "Trend". しっかし、すごい引数名でびびるわ... nullになにかを代入する日が来るとは...
- lshort: 値は論理値。TRUEのとき、truncation lagパラメータが trunc(4+(n/100)^0.25), FALSEのとき trunc(12+(n/100)^0.25) となるのだそうだ。なんだかよくわからん、KPSS検定についてちゃんと勉強しなきゃ。
pp.test(x, alternative, lshort): H0「xは単位根を持つ」のPP検定。定数項と線形トレンドを含む。引数は
- x: 単変量時系列
- alternative: 対立仮説の指定。adf.test()と同じ。
- type: 検定のタイプ。値は"Z(alpha)", "Z(t_alpha)". ううむ、PP検定について不勉強で、なにいってんだかさっぱりわからん。
- lshort: kpss.test()と同じ。
urcaパッケージ.
単位根検定のためのパッケージとして最も有名なのはこれだと思う。
ur.df(y, type, lags, selectlags): ADF検定。引数は
- y: 時系列。
- type: "none"だと切片もトレンドもなし, "drift"だと切片をいれる、"trend"だと切片とトレンドをいれる。これは差分時系列のモデルのことをいってんですよね。
- lags: ラグ次数.
- selectlags. "Fixed"だとなにもしない。"AIC", "BIC"だとlagsに指定した以下の次数をAICなりBICなりで勝手に選んでくれる。
ur.ers(y, type, model, lag.max): Ellot, Rothenberg, & Stock の単位根検定だそうだ。なにそれってかんじだけど、どうやらADF-GLS検定のことらしい。ADF検定より検定力が高いって前に読んだことがあるぞ、理屈はさっぱり理解できなかったけど。引数は
- y: 時系列.
- type: 値は"DF-GLS", "P-test". 後者は誤差項の系列相関を考慮するのだそうだ。へー。
- mode: 値は"constant", "trend". これって元の時系列の話？それとも差分時系列の話？
- lag.max: よくわからんけど、ラグ次数の最大値だそうな。
ur.kpss(y, type, lags, use.lag): KPSS検定. 引数は
- y: 時系列.
- type: 値は"mu", "tau". 前者は切片のみ、後者は切片と線形トレンド。
- lags: 値は"short", "long", "nil"。"short"にするとラグ次数はほにゃらら(メモ省略)、"long"にするとラグ次数はほにゃらら、"nil"にする誤差項の指定なし。
- use.lag: lagsを使わずに、ここで最大ラグ次数を指定してもよい。
ur.pp(x, type, model, lags, use.lag): PP検定。引数は
- x: 時系列.
- type: 値は"Z-alpha", "Z-tau". なんだかわからん。
- model: 値は"constant", "trend".
- lags: 値は"short", "long". マニュアルに説明が書いてない...
- use.lag: lagsを使わずにここで次数を指定してもよい。
ur.sp(y, type, pol.deg, signif): Schmidt & Phillips単位根検定。黒住(2008)ではADF検定ではない「その他の検定」というところで名前だけ出てくる。いわく、ADF検定ってのはWaldタイプの検定なんだけど(なるほど、最尤推定量の差の検定なわけね)、これはラグランジュ乗数検定なんだってさ。そういわれても困るけどな！
引数は
- y: 時系列
- type: 値は"tau", "rho". そういう種類があるんだってさ。しらんがな。
- pol.deg: 多項式の次数だそうだ。値は1,2,3,4。
- signif: 有意水準。値は0.01, 0.05, 0.1.
- y: 時系列
- model: 値は"intercept", "trend", "both". potential breakが切片で起きるか、線形トレンドで起きるか、両方で起きるか。
- lag: ラグ次数の最大値。指定しなくてもいいらしい(えっ、どういうこと？)
forecastパッケージ
泣く子も黙る(?) 有名パッケージ。中の人Hyndman先生は、ただいまこのパッケージのtidyverse対応版であるfableパッケージを鋭意ご製作中らしいのだが、単位根検定関係はまだ移植していない模様。
ndiffs(x, alpha, test, type, max.d, ...): 超お手軽な単位根検定の関数。実はこれ、昨年社外で時系列解析のセミナーやった時にはじめて存在を知り、あまりの簡単さにがっくり膝をついた次第である。だってあれですよ、検定統計量とか一切無視で、単に「何階差分をとるべし」という整数値だけをぽろっと返してくるんですよ。それって、カツカレーくださいって言ったらカツカレーが出てきたようなものじゃないですか...(ちょっとちがうか)
引数は
- x: 時系列
- alpha: 有意水準. 0.01から0.1までの値。
- test: 値は"kpss", "adf", "pp"。
- type: 値は"level", "trend". 恥ずかしながら、これがよく理解できてないんですが... 差分時系列に切片だけはいるか(つまりドリフトつきモデルか)、$t$の1次項もはいるか(つまりトレンドつきモデルか)、ということでしょうか。ってつぶやいてないで、今度ちゃんと調べよう。
- max.d: ラグ次数の最大値。AICかなんかで自動選択してるんだろうなあ。これも今度調べてみよう。
- ...: 単位根検定に渡す引数。たぶんurcaパッケージに渡るんだと思う。
[02/16 追記: 上記の想像通りで、ndiffs()のコードをざっと眺めたところ、ur.df(), ur.kpss(), ur.pp()のいずれかをコールしている模様。次のような対応になっているようだ：
ndiffs(x, test="adf",type="level") → ur.df(x, type="drift")をコール
ndiffs(x, test="adf",type="trend") → ur.df(x, type="trend")をコール
ndiffs(x, test="kpss", type="level") → ur.kpss(x, type = "mu")をコール
ndiffs(x, test="kpss", type="trend") → ur.kpss(x, type="tau")をコール
ndiffs(x, test="pp", type="level") → ur.pp(x, type="Z-tau", model="constant")をコール
ndiffs(x, test="pp", type="trend")→ ur.pp(x, type="Z-tau", model="trend")をコール]
nsdiffs(x, alpha, test, max.D, ...): ndiffs()のSARIMA版、つまり、季節について何階差分をとるべきかを返す。testの値は"seas", "ch", "hegy", "oscb"。おっと、hegyってのがある... uroot::hegy()と同じことなのかなあ...
そのほかのパッケージ
uroot::hegy(x, ...): Hylleberg, Engle, Granger & Yoo の季節単位根検定、だそうだ。勉強不足であれですけど、要するにモデルに季節ダミーが入りますってこと？ それとも、SARIMA(p,d,q,P,D,Q)モデルを当てはめるときにDをどうするかってこと？ いやまてよ、それって実は同じことなの？ うーん、よくわからん。引数はいっぱいあるので省略する。
CADFtest::CADFtest(model, x, type, data, max.lag, min.lag.y, min.lag.X, max.lag.X, dname, criterion, ...): Covariate ADF単位根検定、だそうだ。どうやら、共変量をいれたADF検定らしい。へええ！そんなのあるんだ！ ってことはあれですか、時系列変数間で回帰するとき、yが単位根を持っていてもこの検定に通るなら差分取らなくていいよってことっすか？ それ助かるわー... J.Stat.Softwareに論文が出てるらしい、今度読んでみよう。
MultipleBubbles::ADF_FL(y, adflag, mflag): ラグ次数を固定したADF検定。adflagがラグ次数、mflagは1のとき切片のみ、2のとき切片とトレンド、3のとき両方なし。これ、関数のヘルプを見る限りふつうのADF検定なんだけど、このパッケージは「バブルの存在をチェックするためのADF検定」のパッケージなので、きっとなんか特別なことをやっているんだろうなあ。
MultipleBubbles::ADF_IC(y, adflag, mflag, IC): 上と同じだが、mflagは最大ラグ次数となり、IC=1とするとAIC, 2とするとBICで次数選択するらしい。
fUnitRootsパッケージ: チューリッヒのThe Rmetrics Accoc.というところが本やRパッケージをたくさん出していて、これもそのひとつ。ADF検定の関数を自前で持っている(unitrootTest(), adfTest())。ほかにurcaパッケージへのラッパーもある。
FinTS::Unitroot(x, trend, medhod, lags): これはTsay(2005)という本のコンパニオン・パッケージ。この関数はfUnitRootsパッケージへのラッパーだそうだ。
aTSAパッケージ: SASとよく似た出力を返す時系列分析パッケージなのだそうだ(ははは)。adf.test(), pp.test(), kpss.test()を持っている。
Rの時系列分析パッケージとしては、ほかにTSAという強力なやつがあるんだけど、単位根検定の関数は持っていないようだ。