elsur.jpn.org >

2020年4月 1日 (水)

実にどうでもいい話ですが... すいませんこれは私による私のためのメモです。

年明けに、音楽学者・磯山雅の「マタイ受難曲」という本を読んだのがきっかけで、妙にバッハのマタイ受難曲にはまってしまい、仕事しながらずーっと聴いていた。ついには柄にもなく、たまたまみつけたコンサートのチケットまで取ってしまったのだが(クラシックのコンサートなんてどれだけプチブルなのかとびびりつつ)、そこでコロナ禍が来襲し、コンサートはあえなく中止、のみならず自宅から出られないという事態に陥った。

残念だなあ、とぼんやりyoutubeを眺めていたら、全編演奏の動画が結構あるので驚いた。ひとが口あけて歌っているのをみてなにが楽しいのかと思いますが、これが妙に楽しいんですね。冴えない感じのおっさんが天使のような歌声の持ち主だったりして。

というわけで、youtube上でみつかった、バッハ「マタイ受難曲」の全編演奏の動画を挙げておく(音楽のみで画像は静止画、というのは除外)。正確な演奏年はわからないのが多いんだけど、だいたい新しそうなのから順に。

いきなり再生回数の少ないやつから始まるけど、マレーシア・バッハ・フェスティバル・オーケストラ、David Chin指揮。2019年の演奏らしい。英語と中国語字幕付きというちょっと珍しいライブ動画である。1曲目の合唱の出だし、日本語で「来なさい娘たち、ともに嘆きましょう」は「来吧、女兒們、来興我們哀悼」なのだそうです。

オランダバッハ協会、指揮はヨス・ファン・フェルトホーフェン、福音書記者はベンジャミン・ヒューレットというイケメン。2019年投稿だから、2018-2019年ごろのライブであろう。映像としていちばん見ごたえがあるのはこの動画だと思う。ネットでの発信に力をいれているようで、他にもいろいろ説明が充実してそうだ。

イングリッシュ・バロック・ソロイスツ, 指揮はジョン・エリオット・ガーディナー。おそらく2016年ごろの演奏。福音書記者はマーク・パドモアという、他の動画にも出現するひげのおっさんで、音楽の素養もドイツ語の知識も全くない私だが、このおっさんの謎の説得力にはびびってしまう。この人が新聞の勧誘に来たら取っちゃうかも。

コレギウム・ヴォカーレ・ゲント, 指揮はフィリップ・ヘレヴェッヘ, 福音書記者はクリストフ・プレガルディエンというおっさん。演奏年は不明だが2018年投稿。英語の字幕付き。

Hofkapelle Munchenってかいてあるんだけど、ホーフカプレ?ってどういう意味なの? 指揮はChristian Fliegnerって人、福音書記者はBenjamin Glaubitzというお兄さん。2017年投稿。

Verbier Festival Chamber Orchestra, Thomas Quasthoff指揮, 2015年。ヴェルビエってのはスイスのリゾート地らしい。指揮者が椅子に座ってるんでどうしたのかと思ったら、この方はもともとサリドマイド薬害で障害をもっていて、かつては有名な声楽家だったんだけど引退したんだそうな。福音書記者は... おっと、またマーク・パドモアだ。

説明がロシア語なのでさっぱりわかんないんだけど、機械翻訳によれば、「バッハ・アンサンブル」の演奏で、指揮はヘルムート・リリング。福音書記者のメガネのおっさんはローター・オディニウスという人かしらん? 2014年。

ロイヤル・コンセルトヘボウ管弦楽団, 指揮はイヴァン・フィッシャー、福音書記者はまたもマーク・パドモアさん、2012年。DVDやBlu-rayにもなっているらしい。このメモをとるために各動画の1曲目だけをそれぞれぼーっと聞いてたんだけど、この演奏にはすっかり引き込まれてしまい、2曲目、3曲目と進んで途中で止めようにもなかなか止められない、という感じであった。

アムステルダム・バロック管弦楽団, トン・コープマン指揮。2005年、2日に分けてやった模様。英語字幕付き。

ライプツィヒ・ゲヴァントハウス管弦楽団、指揮はゲオルク・クリストフ・ビラー、福音書記者はマルティン・ペッツォルトという人、1998年。音質がちょっと悪いかな?

サイトウ・キネン・オーケストラ、小澤征爾指揮, 福音書記者ジョン・マーク・エインズリー、1997年。ありがたいことに日本語字幕付きである。

ブランデンブルク・コンソート, スティーヴン・クレオバリー指揮, 福音書記者はRogers Covey-Crumpという人。1994年。

ミュンヘン・バッハ・コレギウム, エノッホ・ツー・グッテンベルク指揮。福音書記者はクレス=ホーカン・アーンショーという人。2015年投稿だがかなり古い演奏だと思う。

よくわかんないけど、指揮はニコラウス・アーノンクールらしい。1985年。残念ながら音質も画質もよくない。

ミュンヘン・バッハ管弦楽団、カール・リヒター指揮、1971年。すごい、こんな映像があるんですね。英語字幕つき。調べてみたら、どうやらDVDになっている模様

というわけで、この危機的事態がいつまで続くのかわからないが、バッハの長い曲など聞きつつ、なすべきことをなし、事態が収束するか自分も感染して死ぬのを待つことにしよう。

雑記 - youtubeでみられる「マタイ受難曲」

2020年3月22日 (日)

 Tポイントとかdポイントとか、ああいうロイヤリティ・プログラムの「ポイント」についての資料。仕事の関係で目を通した。論文じゃないけど、読んだものは何でも記録しておこう、ということで...

翁百合(2019) ポイント経済化について -マクロ経済や金融システムへのインプリケーションを探る-. 日本総研リサーチレポート, 2019-010.

 概観:

  • NRIによれば2018年度のポイント発行額は9546億円。これは企業の売上とかにポイント適用率と還元率をかけて推計している。いっぽう金融庁は引当金(負債に属する)の合計額の推移を調べている。
  • クレカ・家電量販・携帯電話で発行額が多い。伸びているのはプラットフォーマー企業のポイント。
  • ポイントの利用率(失効率)はよくわからない。
  • ポイント交換のネットワークは複雑化している。必ずしも交換が進んでいるわけではない。交換自体をビジネス化する動きもある(ネットマイルとか)。
  • プラットフォーマーはキャッシュレス決済の経済圏をつくりポイント付加で囲い込みを図っている。
  • ポイントを投資に転換する動きもある。これについては金融庁が前向きな方針を示している [←へー]

 ミクロにみると、個社発行ポイント(プラットフォーマーでない企業のポイント)とは、企業が消費者の価格弾力性に注目し、将来の値引き分の決済手段をトークンとして渡すもの。いっぽうプラットフォーマーは、還元額の多くを加盟店に負担させる形でトークンを付与する、というちがいがある。

 プラットフォーマーがポイント還元を積極的に活用するとなにがおきるか? 3つにわけてみていく。

 A. 消費の質の変化。先行研究ではふたつの視点がある。
 A1. ポイントがスイッチングコストとして作用し、ロックインが起きる。この観点からいえばポイント交換はスイッチングコストを下げる。
 A2. 行動経済学的視点でいうと、ポイント付与の表現を工夫したり参照点を変えたりしてフレーミング効果を起こす。[このくだりは先行研究列挙に近い。寺地・近(2011行動経済学), 泉谷(2013放送大の紀要), 中川(2015行動経済学), 山本(2011立教の紀要), Liu(2007 J.Mktg)]
 というわけで、ポイントはロックインにより消費者余剰を低下させている可能性がある。また、キャッシュレス化、データ活用による個別消費提案、消費者行動の非合理化(ポイントを貯めるための貯めるための追加消費とか)が起こりうる。

 B. 価格のカスタム化。先行研究ではポイントを非線形プライシングととらえる見方がある。また、企業は景気が悪くなると囲い込んだ顧客に高めの価格をつける傾向があるので、counter cyclicalな効果があるかもといわれている。
 というわけで、物価指数が実態から離れるかも。なお、プラットフォーマー型ポイントは価格決定力が加盟店にないからcounter cyclicalな効果は小さかろう。

 C. ポイントの疑似貨幣化。先行研究では...[中略]。ただし、仮に年5000億のポイントが決済に使われているとしても、電子マネーの決済金額は約5.5兆、それほど大きくない。また、1単位当たりの価値がばらばらだし交換が限定されているし期限があるし還流しない。貨幣の3機能(価値尺度、流通手段、価値貯蔵)をすべて満たさない。今後のプラットフォーマー型ポイントの進展によってはわかんないけど。
 むしろポイントを貯めること自体に楽しみを見出す人も多いわけで、法定通貨にない魅力があるのかもしれない。貨幣の多様性の展望を拓くものともいえる。
 云々。

 ... いっちゃなんだけど、もう少し体裁に気を配ればいいのに、と思った。ちゃんと校正するとか、数式を数式フォントで書くとか。こういうのを読む専門家の方は、そのへんあんまし気にしないのかなあ...

 ポイントについての研究を探すとみつかるのはたいてい購買をアウトカムにとった行動経済学的な話なんだけど(ポイント付与と値引きがどう違うかとか)、そうじゃなくて、ポイント獲得行動そのものといえばいいのか、ケインズのいう貨幣愛ならぬポイント愛といえばいいのか、この資料でいうとCの最後の論点に関心があって、理論面での先行研究を探しているんだけど、見当たらない。このレポートでもアンケートを紹介しているだけだった。ううむ。

論文:マーケティング - 読了:翁(2019) マクロ経済に対する「ポイント」の影響

2020年3月19日 (木)

 コロナ騒動で気が休まらない毎日でありますが、感染者数の推移や予測がメディアにさかんに取り上げられていて、ああ、こういう研究ができたら面白かったろうな... などと無責任なことをぼんやり考えていた次第である。
 で、フォルダを整理していたら、「いずれ読む」PDFのフォルダにそっち方面の概説論文がはいっているのをみつけた。全然覚えていないんだけど、仕事の参考になるかと思って深く考えずに保存していたのだと思う。
 せっかくなので、仕事の待ち時間にめくってみた。

西浦博, 稲葉寿 (2006) 感染症流行の予測:感染症数理モデルにおける定量的課題. 統計数理, 54(2), 461-480.

 まずはSIRモデルから。
 感染する可能性のある人を$S$, 感染している人を$I$, 免疫を獲得したか死んだ人を$R$として、
 $\frac{dS(t)}{dt} = -\beta S(t) I(t)$
 $\frac{dI(t)}{dt} = \beta S(t) I(t) - \gamma I(t)$
 $\frac{dR(t)}{dt} = \gamma I(t)$
$\beta$は感染率。$\beta I(t)$を感染力という。$\gamma$は回復率や隔離率で、$\gamma^{-1}$を感染期間が指数分布に従うと仮定した平均値とする。短期的な流行のモデルなので人口動態は無視する。

 総人口を1として、平衡状態$(S(t), I(t), R(t)) = (1,0,0)$において
 $\frac{dI(t)}{dt} = (\beta - \gamma) I(t)$
が成り立つ。感染初期に感染者は$I(t) \approx I(0) \exp{(\beta -\gamma)t}$と指数関数的に増えるわけだ。病気が集団に侵入可能になるのは$\beta - \gamma > 0$のときである。ここから
 $R_0 = \frac{\beta}{\gamma}$
を基本再生産数という。R noughtと発音する。[←完璧な門外漢なので、このくだりを読んだだけでもう元をとったという感じだ。SIRモデルについては読んだことがあったけど、いまメディアでみかける基本再生産数って全然意味がわかってなかった]

 時刻 $t$までに感染を経験した人を $N(t) = S(0) - S(t)$とすると... [残念ながらこの段落は知識不足で理解できなかったので、メモも省略。シグマ代数ってなに?]

 SIRモデルの拡張。
 $I$をさらに、曝露したけどまだ感染性を持たない状態$E$と感染性を持つ状態$I$にわけると
 $\frac{dS(t)}{dt} = -\beta S(t) I(t)$
 $\frac{dE(t)}{dt} = \beta S(t) I(t) - \epsilon E(t)$
 $\frac{dI(t)}{dt} = \epsilon E(t) - \gamma I(t)$
 $\frac{dR(t)}{dt} = \gamma I(t)$
$\epsilon$は「曝露後に感染性を得る率」。これをSEIRモデルという。
 ほかにもモデルが作れる。免疫が獲得できないSIRSモデルとか。

 基本再生産数の意義として、予防接種率の達成目標について考えると...
 効果$\epsilon$のワクチンを接種率$p$で実施したとき、免疫を持たない人の割合は$1-\epsilon p $だから、接触パターンがhomogeneousなら再生産数は$(1-\epsilon p) R_0$。これを1未満にすれば感染症は終息する。$p$について解くと $p > 1 / \epsilon (1-1/R_0)$だ。[←うぉー おもしれええ]
 $R_0$の推定は、さすがに決定論的に考えてると精度とかを考慮できないので、確率論的なモデルを考えて...[めんどくさいからメモは省略するけど、そんなにものすごく難しい話ではなさそう。ふつうはML推定量が閉形式で書き下ろせるらしい。$\beta$の時間変動とか感受性の異質性とかを考慮したモデルをMCMC推定したりするのもあるのだそうだ。ふへえええええ]

 実際には伝播能力は行動の変化とか公衆衛生施策とかで時間変動するので、そういうのを考慮したのを効果的再生産数($R_t$)という。実際の感染ネットワークから求めたり、理論的に求めたりできる。SARSの流行時にWHOが緊急事態のアラートを出した直後に$R_t$が1未満になったという研究があるそうだ[←へえええええ!]

 感染待ち時間・感染性期間の分布を天下りに与えるのではなく、観察に基づき適切な分布を推定するというアプローチもある。
 感染待ち時間は潜伏期間と似ているので(潜伏期間の途中から感染性があることが多いから、ほんとはちょっとちがうんだけど)、調査で調べた潜伏期間の分布で代用したりする(対数正規分布に従うことが多い)。検疫に必要な目安の時間がわかったりする。[途中から難しくなってきたので大幅中略]
 いっぽう感染性期間の分布を知るのは難しい。感染待ち期間と潜伏期間は似ていることが多いけど、感染性期間と症候性期間はちがう。たとえばHIVの場合、感染した直後にウィルスが多く、いったん下がってAIDS発症の前にまた上がる。[...] 実際のウィルス量と疫学的な感染性との評価指標との相関はまだ確実でなく、今後の課題である。[←へえー。個人レベルの機序と集団レベルの現象の関係が意外にはっきりしないってことなのね... ウィルス量の増減も人によっていろいろってことかしらん]

 接触性の異質性を考慮するという方向の拡張もある。[...中略するけど、ひとつのアプローチはそのなかでは接触パターンが均一であるような複数のsubpopulationを考えるという方法で、これをmultitype epidemicモデルというのだそうだ。よくわかんないけど潜在クラスモデルみたいになるのだろうか...]
 最近関心が集まっているのは、接触関数の頻度分布がべき乗分布に従うとき(スケールフリーなとき)で、感染症の制圧は困難を極める。
 接触頻度に関する定量的説明は性感染症から明らかになるだろうと期待されている。接触回数の計数が観察可能だから。[←これ、誰と何回セックスしたかインタビューで正直に答えるだろうから、って意味だろうか...?]

 あれこれ考えているとどんどん複雑になるので、観察データじゃなくて介入研究でパラメータ推定することも多い。[←うぉー面白い!と思ったんだけど、どういう実験になるのか想像がつかなかった...]

 云々、云々。

 ...半分くらいはちんぷんかんぷんだったんだけど、仕事と直接関係しないという気楽さもあって、いやー、面白かった。もっと頭がよくて若かったら、こういうのやりたかったなあ... (← 一時の気の迷い)
 ともあれ、研究者の方々のご奮闘をお祈りいたします。

論文:データ解析(2018-) - 読了:西浦・稲葉(2006) 感染症流行の予測

2020年3月11日 (水)

 仕事の都合であれこれ考えていたらわけがわからなくなってしまったので、ちょっとシミュレーションをやってみたら、さらに新たな疑問が生じ... という話を記録しておく。

背景
 ある量的変数の時系列 $Y_t$ があって、それに影響を及ぼしているであろうなんらかの変数の時系列$X_t$があるとする。なんでもいいんだけど、たとえば売上と広告出稿量とか。
 で、$Y_t$に対する$X_t$の効果の大きさをデータから推定したい。こういうこと、よくありますよね。

 話を簡単にするために、分析者は次のことを知っているとしよう。

  • $X_t$の$Y_t$に対する効果は即時的、かつ線形である。つまり、
     $Y_t = \alpha + \beta X_t + V_t$
    とモデル化できる。
  • $X_t$は定常である。なんでそんなことを知っているのかわからんが、なにか実質的な知識があるのでしょう。
  • $Y_t$の撹乱項$V_t$は一次自己回帰過程 (AR(1)過程)に従う。つまり、
     $V_t = \rho V_{t-1} + E_t, \ \ E_t \sim N(0, \sigma_E^2)$
    とモデル化できる。これも、まあ、ありそうな仮定ですわね。時系列の自己相関関数を観察した結果、これAR(1)でいいんじゃね? と思うことはさほど珍しくないと思う。
  • 自己回帰係数$\rho$は$0 \leq \rho \leq 1$である。つまり、$V_t$は非負の自己回帰係数を持つ定常AR(1)過程に従う($0 \leq \rho <1$)か、単位根AR(1)過程に従う($\rho = 1$)かのどちらかである。これもそれほど変な仮定ではないと思う。自己回帰はふつう正だし、もし負だったらデータの観察でそれと気がつくだろう。仮に$\rho >1$なら$V_t$はどっかに飛んでいってしまうはずなので、これも除外できる。でも、$\rho <1$か$\rho = 1$かは簡単には判断できない。

問題
 さて、分析者が関心を持っているのは$\beta$である。ここで疑問なのは、$\beta$の推定誤差は$\rho$とどういう風に関連しているのか、という点である。特に次の点について知りたい。

  • Q1. $\rho$が大きいとき、つまり結果変数の撹乱項が高い自己相関を持つとき、$\beta$の推定誤差は大きくなるか、それとも小さくなるか。
  • Q2. $\rho = 1$であるのにそれと知らず、$\rho < 1$と仮定するモデルを推定したとき、$\beta$の推定誤差は大きくなるか。

 いずれも、私にとってはちょっと切実な疑問だ。
 Q1についていえば... 仕事のなかで時系列データの分析が生じる際、まずは目的変数の時系列を観察して、これからわざわざ説明変数のデータを集めモデリングする労力は報われるかしらん? と算段することが多いと思う。そういう場面での手掛かりがほしい。見通しが暗いなら早めに白旗を上げたい。
 Q2についていえば... 理屈の上からいえば、$\rho=1$かそうでないかでAR(1)過程の挙動はがらっとかわる。でも、実データの背後にあるデータ生成プロセスが$\rho=1$かどうかなんて、どうせわかりっこない。だから、きっと$\rho<1$だよねと信じてAR(1)誤差を仮定したり、いや$\rho=1$にちがいないと信じて差分時系列を分析したりするわけである(そんなことないですか?)。でも心には常に不安の影が付きまとう。$\rho$についての仮定をしくじることが、$\beta$の推定にとってどのくらい致命的なのかを知りたい。

方法
 気になって夜も眠れないので(大げさな表現)、簡単なシミュレーションをやってみました。

 $t=-100,\ldots,100$について、データを次のように生成した。
 $x'_t \sim unif(0, 1)$
 $e_t \sim N(0, 1)$
 $v_t = \rho v_{t-1} + e_t$
 $y'_t = x't + v_t$
つまり$\beta = 1$である。で、$x'_t, y'_t$から$t=1$以降を切り出し、それぞれを中心化する。これを$x_t, y_t$とする。
 $\rho$を$0, 0.1, 0.2, \ldots, 1.0$の11通りに動かし、それぞれについて1000個のデータセットを生成した。
 それぞれについてデータセットをひとつ選び$x_t, y_t$を描画すると、こんな感じである。

plotData.1.jpg

 モデルの推定方法として、次の6つを試す。

  • A. 単にOLS推定
  • B. Rのforecast::Arima()でML推定
  • C. Rのnlme::gls()でFGLS推定
  • D. 状態空間モデルとして定式化し、RのKFASパッケージを使ってカルマンフィルタで推定
  • E. Mplusでベイズ推定
  • F. Stanでベイズ推定

 推定にあたっては、次の5つのアプローチを試してみる。

  • アプローチ1. $\rho = 0$と仮定して、
     $y_t = \beta x_t + e_t$
     $e_t \sim N(0, \sigma_e^2)$
  • アプローチ2. $|\rho| < 1$と仮定して、
     $y_t = \beta x_t + v_t$
     $v_t = \rho v_{t-1} + e_t, \ \ |\rho| < 1$
     $e_t \sim N(0, \sigma_e^2)$
  • アプローチ3. $\rho = 1$と仮定し、$\Delta y_t = y_t - y_{t-1}, \Delta x_t = x_t - x_{t-1}$について、
     $\Delta y_t = \beta \Delta x_t + e_t$
     $e_t \sim N(0, \sigma_e^2)$
  • アプローチ4. $y_t$について単位根検定を行い、単位根を持たないと判断されたらアプローチ2, 持つと判断されたらアプローチ3に進む。
  • アプローチ5: アプローチ2と同じなんだけど、$\rho$の範囲について仮定しない。

 私の乏しい理解では、A(OLS), B(ML), C(FGLS)では撹乱項の定常性が仮定されるので、アプローチ5は選べない。いっぽうD(カルマンフィルタ), E,F(ベイズ)では、$\rho$を明示的に制約しないかぎりアプローチ5になる... というように理解しているのだけれど、ここ、全然自信がない。悲しい。誰か私に教えてくださらないでしょうか。

 まあいいや!とにかくシミュレーションしてみたのであります。

単位根検定
 まず、アプローチ4で必要になる単位根検定の結果について紹介しておく。Dickey-Fuller検定を使った。Rコードはこんな感じ。$y_t$がdfIn$y_centeredにはいっている。

library(urca)
oDFTest <- ur.df(dfIn$y_centered, type = "none", lags = 1)
nD <- if_else(oDFTest@teststat[1,1] < oDFTest@cval[2], 0, 1)

トレンドもドリフトもないAR(1)であると知っているので、type ="none", lags=1と決め打ちしている。$H_0: \rho=1$が 5%水準で棄却されたら$\rho<1$と判断し、そうでなかったら$\rho=1$と判断することにして、後者の場合にnDを1としている。

 素朴に考えると、本当は$\rho=1$であるデータセットのうち95%くらいが$\rho=1$と判断されてほしい($H_0$が真のときに誤って棄却される確率は5%であってほしい)。本当は$\rho<1$である場合は、なるべく多くのデータセットが$\rho<1$と判断されてほしい。
 さて、11通りの$\rho$の、それぞれ1000個のデータセットのうち、$\rho = 1$と判断されたデータセットの数は?

plotDFtest.1.jpg

 蓋をあけてみると... $\rho=1$であるデータセットのうち、$H_0$が棄却されなかった($|\rho|=1$と判断された)データセットは95%には到底及ばず、実に68%にとどまる。なぜだろう? よくわからない。
 $\rho<1$であるにもかかわらず$H_0$が棄却されないデータセットは、当然ながら$\rho$が1に近づくにつれて増えるのだが、それでも$\rho = 0.9$のときに9%。$\rho = 0.8$までではほとんど生じない。
 このように、DF検定の場合、保守側($H_0: \rho = 1$が棄却されない側)に大きくバイアスがかかるようだ。$H_0$を反対側に設定するPP検定であればまた違う結果になるだろうけど...
 まあとにかく、ここで確認しておきたいのは、$\rho = 1$かそうでないかなんて、そんなに簡単にはわからないよね、ということである。

選手紹介
 お待たせしました、選手入場です。拍手でお迎え下さい。

 A. 単にOLS推定。アプローチ1($\rho=0$と仮定)のRコードは

oModel <- lm(y_centered ~ 0 + x_centered, data = dfIn)

自己相関を華麗にスルーしちゃうわけだが、これ、それほど捨てたもんじゃないと思う次第である。だって、$|\rho| < 1$であれば、$\hat{\beta}$は少なくとも一致推定量ではあるわけでしょう?
 アプローチ1, 3($\rho=1$と仮定)を試してみる。

 B. forecast::Arima()で最尤推定。Rによる時系列モデリングの定番 forecast パッケージはArima()という関数をご用意している(実はstat::arima()へのラッパーである)。アプローチ2($|\rho|<1$と仮定)のRコードはこんな感じ。

library(forecast)
oModel <- Arima(
y = dfIn$y_centered,
include.mean = FALSE,
order = c(1, 0, 0),
xreg = dfIn$x_centered ,
method = "ML"
)

アプローチ3($\rho=1$と仮定)ならorder = c(0,1,0)となる。アプローチ2, 3, 4(DF検定で切り替え)を試す。

 C. nlme::gls()でFGLS推定 ついでにnlme::gls()でFGLS推定を試してみる。計量経済学の教科書に載っているのは、最尤推定じゃなくてこっちのほうですね。
 アプローチ2($|\rho|<1$と仮定)のRコードはこんな感じ。

library(nlme)
oModel <- gls(
y_centered ~ 0 + x_centered,
corr = corARMA(p=1, q=0),
data = dfIn
)

アプローチ2のみ試す。

 D. 状態空間モデルとして定式化し、RのKFASパッケージを使ってカルマンフィルタで推定。おさらいすると、アプローチ2,5のモデルは
 $y_t = \beta x_t + v_t$
 $v_t = \rho v_{t-1} + e_t, \ \ e_t \sim N(0, \sigma_e^2)$
である。これを状態空間表現に書き換えよう。状態変数は$\beta$と$v_t$だと考え、縦に積んで$\alpha_t = [\beta, \ v_t]'$としよう。
観察方程式は
 $Y_t = Z \alpha_t$
ただし $Z = [x_t, 1]$である。撹乱項がないことに注意。状態方程式は
 $\alpha_t = T \alpha_{t-1} + R \eta, \ \ \eta \sim MVN(0, Q)$
遷移行列$T$は2x2の対角行列で、対角要素は$1, \rho$である。$R$は2x2の単位行列で、状態撹乱項は$\eta = [0, e_t]'$とする。共分散行列$Q$は2x2で、右下に$\sigma_e^2$がはいり、残りが0。
 こいつをカルマンフィルタで推定する。Rコードはこんな感じ。遷移行列$T$に未知パラメータが入っているので、ちょっと面倒くさい。

library(KFAS)
mgZ <- array(dim = c(1, 2, nrow(dfIn)))
mgZ[1,1,] <- as.vector(dfIn$x_centered)
mgZ[1,2,] <- 1
oModel <- SSModel(
dfIn$y_centered ~ -1 + SSMcustom(
Z = mgZ,
T = matrix(c(1, 0, 0, 1), nrow = 2), # 最後の要素がrho
R = matrix(c(1, 0, 0, 1), nrow = 2),
Q = matrix(c(0, 0, 0, 1), nrow = 2), # 最後の要素がsigma_e
P1 = matrix(c(0,0,0,0), nrow = 2),
P1inf = matrix(c(1,0,0,1), nrow = 2)
),
H = matrix(0)
)
sub_update <- function(par, model){
# par: (sigma_eの対数, rho)
# model: 現在のモデル
model$T[2,2,1] <- par[2]
model$Q[2,2,1] <- exp(par[1])
return(model)
}
oFitted <- fitSSM(
oModel,
inits = c(log(var(dfIn$x_centered)/2), 0.5) ,
updatefn = sub_update,
lower = c(-10, -2),
upper = c(+10, +2),
method = "L-BFGS-B"
)
oEstimated <- KFS(oFitted$model)

計算の都合上$|\rho| < 2$と制約しているものの、このモデルは$|\rho| = 1$という仮定を置かないモデル、つまりアプローチ5ということになると思うのですが... 正しいでしょうか?

 E. Mplusでベイズ推定。ここまではRなどという古くさい言語を使っておりましたが、いけてる分析者なら、ここは当然 Mplus ですよね! (すいません冗談です)
 構造方程式モデリングの世界ではもはや標準となっているソフトウェア Mplus だが、実はN=1の時系列分析についても便利な機能を持っているのである。
 Mplusのコードはこんな感じ。

DATA: 
FILE = est2.dat;
VARIABLE:
NAMES = y x;
MISSING=.;
ANALYSIS:
ESTIMATOR = BAYES;
BITERATIONS = (2000);
MODEL:
v by (&1);
v;
v on v&1;
y on v@1 x;
[y@0];
y@0;

 MODELコマンドがわかりにくいが、1行目は「vは潜在変数ですが指標を持っていません。モデルのなかでラグ1を使わせて下さい」。2行目は「その残差分散(つまり$\sigma_e$)を自由推定したいです」。3行目で$v_t$の自己回帰を定義し(v&1とは$v_{t-1}$を表す)、4行目で$y_t$のモデルを定義している($v_t$の回帰係数は1に固定している)。ほっとくと$y_t$の切片と残差分散を推定してしまうので、最後の2行でそれを抑止している。
 このとき$v_t$の自己回帰係数($\rho$)は範囲が制約されていない、つまりアプローチ5だ、というのが私の理解なのだが、正しいだろうか...?

 F. Stanでベイズ推定。ほんとはここまでやるつもりはなかったんだけど、毒も食らわば皿まで、ということで... Stanファイルはこんな感じ。

data {
int T;
vector[T] y;
vector[T] x;
}
parameters {
real beta;
real<lower=-1, upper =1> rho;
real<lower=0> sigma;
}
model {
vector[T] v;
v = y - beta * x;
v[1] ~ normal(0, sigma/sqrt(1-rho^2));
v[2:T] ~ normal(rho * v[1:(T-1)], sigma);
}

 下から3行目、v[1]でなにか変なことを書いているが、これは$v_{t-1}$が未知であるときの$v_t$の周辺分布の分散が、$v_{t-1}$の下での$v_t$の条件つき分布の分散 $\sigma_e$より大きくなるからである。
 どのくらい大きくなるかというと、
 $v_t = \rho v_{t-1} + e_t$
の両辺の分散をとって
 $Var(v_t) = Var(\rho v_{t-1} + e_t)$
$e_t$と$v_{t-1}$との共分散は0なので、
 $Var(v_t) = \rho^2 Var(v_t) + \sigma_e^2$
ここから
 $Var(v_t) = \sigma_e^2 / (1-\rho^2)$
である。
 というわけで、尤度計算に当たってデータ点ひとつも無駄にしまいという質実剛健な精神に則り、v[1]についても誠意を込めて分布を書いてみたんだけど、でもこれって、$1-\rho^2 \neq 0$、つまり$|\rho| \neq 1$という仮定を暗黙のうちに含んでいないだろうか? そんならいっそ、というわけで、parametersブロックではreal <lower=-1, upper =1> rho; と書いた次第である。
 このモデルとともに、parametersブロックでrhoの範囲制約をなくし、かつ下から3行目を消したモデルも推定してみた。前者はアプローチ2, 後者はアプローチ5に相当すると思うんだけど... うーーん、こういう理解で正しいのだろうか? からきし自信がない。

 まとめると、出場選手は以下の10名である。

  • A1: OLS推定, $\rho = 0$と仮定
  • A3: OLS推定, $\rho = 1$と仮定
  • B2: 最尤推定, $|\rho| < 1$と仮定
  • B3: 最尤推定, $\rho = 1$と仮定
  • B4: 最尤推定, 単位根検定に基づきモデル選択
  • C2: FGLS推定, $|\rho| < 1$と仮定
  • D5: カルマンフィルタ
  • E5: ベイズ推定(Mplus)
  • F2: ベイズ推定(Stan), $|\rho| < 1$と仮定
  • F5: ベイズ推定(Stan)

推定の様子
 この10人の選手に、11水準の$\rho$の各1000個のデータセットについて、$\beta$を推定させた。ただしStanは時間がかかるので、1000個のうち100個だけについて推定するだけで勘弁してやった。

 いったいなにをやっているのか、自分でもよくわかんなくなってきたので、選手B2(最尤推定, $|\rho| < 1$と仮定)による推定結果を図にしてみた。

plotB2.2.jpg

 点はデータセット、横軸は$\beta$の推定値, 縦軸は$\rho$の推定値である。真値が赤で表現してある。こうしてみると、真の$\rho$が1に近づくにつれ、$\rho$の推定値は0に向かって歪むみたいですね。横軸と縦軸の間にはあまり相関がなさそうだ。
 さて、関心があるのは、$\beta$の真値すなわち 1 と、$\beta$の推定値$\hat{\beta}$とのずれである。その大きさを次の2つの指標で評価しよう。

  • バイアス。$\hat{\beta}$の平均から1を引いた値。0に近くないと困る。
  • RMSE、すなわち$(\hat{\beta}-1)^2$の平均の平方根。0に近いほうがありがたい。

結果
 大変ながらくお待たせいたしました!お待たせしすぎたかもしれません! (ちょっと懐かしい冗談だ)
 結果発表です!

plotEstStat.1.jpg

 選手F2, F5のみ試行回数が100である点に注意。
 シンプルなチャートだが、なかなか情報量が多いので、順にみていこう。

 まずは、$\rho = 0$と勝手に想定し、OLS回帰をやってしまった場合。
 真の$\rho$がそれほど大きくなければバイアスはそれほど大きくない。しかしRMSEは$\rho$とともに増大する。$\rho$が1に近づくと、バイアス・RMSEともに急増する。なるほど、撹乱項が非定常に近づいているわけだから、パラメータがうまく推定できないのも道理である。
 やっぱりあれですね、時系列を分析する際には、きちんと自己相関を考慮することが大事ですね。

 次に、$|\rho|<1$と決め打ちした場合(図のパネル2)と、$\rho=1$と決め打ちした場合(パネル3)について。バイアスはそれほど大きくないようだ。RMSEのみ拡大して示す。なお、選手F2はあとで観察することにして、ここでは省略する。

plotEstStat.3.jpg

 選手B2とC2, 選手B3とA3はぴったり重なってしまっている。つまり、$|\rho|<1$と決め打ちする場合、最尤推定するかFGLS推定するかにはたいしたちがいがないわけだ。$\rho=1$と決め打ちする場合に差分をとって最尤推定するかOLS推定するかも... そりゃそうか、きっと推定値は同じだ。
 $\rho$が小さいとき、$\rho=1$と決め打ちするとRMSEが大きくなる。誤ったモデルを指定したわけだから、これは当然である。

 これをみて大変に意外だった点がふたつある。
 第1に、$\rho$が大きくなるにつれ、つまり撹乱項の自己相関が大きくなるにつれ、$\beta$のRMSEが小さくなっているという点である。$e_t$の分散が同じなら$\rho$が大きいほど撹乱項$v_t$の分散は大きくなるんだから、直感的には、$\beta$の推定はより難しくなるんじゃないかと思ったんだけど... なぜだろう???
 もうひとつ、ふへええと言葉にならないため息をついたのは、$\rho$が1に近づくほど、$|\rho|<1$と仮定して定常AR(1)撹乱項を持つ回帰モデルを推定しても、$\rho=1$と仮定して差分時系列について回帰モデルを推定しても、違いがなくなってしまう点である。さきにみたように、単位根検定がしくじりやすいのは$\rho$が1に近いときである。そういうときほど、実は間違えたところでたいした実害はないわけだ。
 このように、単位根検定でモデルを切り替えようが、頭から$|\rho|<1$と信じこんで定常AR(1)撹乱項を持つ回帰モデルを使い続けようが、ほとんどちがいがないようだ。たとえ$\rho=1$だったとしても、である。
 ええええ? わざわざ単位根検定をやったのに、意味なかったわけ...?

 最後に、$\rho$の範囲について制約しないモデル。
 元の図にもどると、選手E5がとんでもないバイアスを持っていることがわかる。どうしたんだMplus! しっかりしろ! 医者だ、医者を呼べ!!
 ... 冗談はともかく、ひょっとすると私がコードを間違えているのかもしれない。残念ながらMplusくんは休場とし、他の選手についてRMSEを拡大してみよう。
plotEstStat.4.jpg

 選手D5 (カルマンフィルタ) について。参考のために選手B2 (最尤推定) と並べてみた。ほとんど変わらない。
 恥かしながらわたくし、撹乱項の定常性を仮定しているB2は$\rho=1$に近づくとだめになるが、定常性を仮定していないD5はうまくいく... という結果になるかな、と思っておりました。不明を恥じる次第であります。へええ、そうなのかー。
 選手F5(Stanでベイズ推定)について。D5との優劣ははっきりしない。参考のためにF2 ($|\rho|<1$と仮定してStanでベイズ推定)と並べてみたところ、$\rho$が小さいところではF2が僅差で勝つが(そりゃそうだ、F2は正しい制約を追加しているわけだから)、$\rho$が1に近づくと僅差で負けるようだ。

まとめ
 というわけで、シミュレーションの結果わかったことをまとめておくと、

  • Q1. 時系列回帰モデル$Y_t = \alpha + \beta X_t + V_t$において、撹乱項$V_t$が高い自己相関を持つとき、$\beta$の推定誤差は大きくなるか、それとも小さくなるか。→小さくなる。
  • Q2. $V_t$の自己回帰係数が1であるのにそれと知らず、1未満と仮定するモデルを推定したとき、$\beta$の推定誤差は大きくなるか。→ならない。

 ううむ。どちらも意外な結果であった。勉強が足りないようだ。

 それにしても、単位根検定の意味ってなんだろう? 時系列分析の教科書には必ず書いてあるけど、やっても意味なくない? ... と思ってしまったのだが、これは私が「時系列回帰で回帰係数を推定する」という場面だけに焦点を絞っているからで、たとえば予測に関心があるなら話はちがうのかもしれない。
 それに、ここでは説明変数時系列が定常であると知っている場合について考えているが、もしそうでなかったら、そりゃあまあ単位根があるかどうか知りたいですわね、みせかけの回帰が怖いから。そういう意味でも、単位根の有無を調べることは、やはり大事なのでありましょう。

 なお、この記事のために書いたコードはすべてGithubにアップしております。自己満足もいいところだがな!

2020/03/16追記: 先週書いたこの記事を見直して、はっと気が付いたので記録しておく。
 お題は次の通りである。
 $y_t = \beta x_t + v_t$
 $v_t = \rho v_{t-1} + e_t, \ \ e_t \sim N(0, \sigma_e^2)$
というAR(1)誤差回帰モデルで、$\rho$が大きくなるほど$\hat{\beta}$のRMSEが小さくなるのはなぜか?

 2本目の式をラグ演算子$L$を使って書き直すと
 $v_t = \rho L v_t + e_t$
 $(1-\rho)L v_t = e_t$
ラグ多項式$(1-\rho)L$を$R(L)$とし、
 $R(L) v_t = e_t$
と書くことにする。$|\rho| < 1$なら$R^{-1}(L)$が定義できて
 $v_t = R^{-1}(L) e_t$
1本目の式に代入して
 $y_t = \beta x_t + R^{-1}(L) e_t$
両辺に$R(L)$をかけて
 $R(L) y_t = \beta R(L) x_t + e_t$
これは($\rho$を既知とすれば)通常のOLS回帰と同じなので、$x_t^* = R(L) x_t$と略記して
 $\displaystyle Var(\hat{\beta}) = \frac{\sigma^2_e}{\sum(x_t^* - \bar{x}^*)^2} $
である。つまり、AR(1)誤差回帰モデルの回帰係数の標準誤差は、ふつうの単回帰のように「撹乱項の分散と独立変数の偏差平方和の比」なのでなく、「AR(1)誤差の裏にあるホワイトノイズの分散と、$x_t - \rho x_{t-1}$の偏差平方和の比」なのである。だから、自己回帰係数が大きいほど分母は大きくなり、$\hat{\beta}$のSEは小さくなるわけだ。うっわー。

plotEstStat.5.jpg

 オレンジ色は、選手B2による$\hat{\beta}$の、真値($\beta=1$)に対するRMSE (再掲)。青色は、各試行について真の$\rho$を既知として$Var(\hat{\beta})$を求め、試行を通じて平均して平方根をとった値。
 な・る・ほ・ど...

雑記:データ解析 - 時系列回帰の撹乱項が自己回帰していると、いったいなにがどうなるのか(追記あり)

2020年2月29日 (土)

 Rをつかっているとformulaを書かないといけないことがあるけど(回帰系の関数を呼ぶときとかに)、なんか便利な書き方がありそうなのに全然使っていない。このたび、ふとstats::formula()のドキュメントを眺めていたら、恥ずかしながら「へぇー」と思うことがいくつかあった。
 この話題、Rを使っている人にとってはたぶん馴染み深い話で、webをちょっと検索するだけで、丁寧に解説しておられるブログ記事をいくつもみつけることができる(たとえばこちら:[R] 予測モデルを作るには formula を活用せよ)。そういうのをきちんと読んできちんと理解しておけばいいんですけどね。そのときは「へぇー」と思うんだけど、すぐに忘れちゃうのである。
 
 というわけで、今日は忘れないようにメモしておこう。なぜか執事風の丁寧語で。お嬢様、羊肉のソテーでございます、って感じで。

  • y ~ model」: 反応 ymodelで指定した予測子でモデル化いたします。
  • a + b」: ab の線形和でございます。
  • a : b」: ab の交互作用でございます。
  • a * b」: abのクロス、すなわち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」: baにネストしております。すなわち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の右辺で01が等価だってのが気持ち悪くてしかたない。責任者でてこい]
  • 「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日 (月)

Bookcover 中国 人口減少の真実 (日経プレミアシリーズ) [a]
村山 宏 / 日本経済新聞出版社 / 2020-01-11
皮肉でもなんでもなく素直に感心しているんだけど、日経新聞の編集委員ともなると、中国メディアの最近の面白記事をピックアップしてコメントをつけただけで、こうして本が一冊書けちゃうんですね。(本当に皮肉ではない。だって誰にでもできることではないじゃないですか)

ノンフィクション(2018-) - 読了:「中国 人口減少の真実」

Bookcover わたしの名は赤〔新訳版〕 (上) (ハヤカワepi文庫) [a]
オルハン パムク / 早川書房 / 2012-01-25
Bookcover わたしの名は赤〔新訳版〕 (下) (ハヤカワepi文庫) [a]
オルハン パムク / 早川書房 / 2012-01-25
年明けから身内の事情で、具体的になにができるわけでもないのにどうも気が休まらない日々が続いていたのだが、新幹線での移動中にずっと手に取っていたのが、以前から気になっていたオルハン・パムクの長編小説であった。良い本は少なからず心の支えとなってくれる。
 いやあ、それにしても、かのノーベル文学賞作家の代表作がこんなエンタメ小説だとは。おそるべしパムク。

Bookcover ポンス・ピラト ほか―カイヨワ幻想物語集 [a]
ロジェ カイヨワ / 景文館書店 / 2013-05-29
私鉄沿線の小さな書店でふと見かけて、R.カイヨワって「遊びと人間」の? 社会学者じゃなかったの? と意外の念に撃たれ、つい手に取った。
 表題作は総督ピラトがイエスを処刑せよと迫られ、友人に相談したりするんだけど(友人はその後の西洋史を丸ごと幻視し、フランスの作家がこのやりとりを小説にするだろうとまでいう)、結局は処刑を拒否、おかげでキリスト教は興りませんでした...という人を食った小説であった。
 ところで、この古風な体裁のペーパーバックを出した景文館書店なる版元は? と疑問に思って検索してみたら、愛知県岡崎市にあるひとり出版社。立ち上げ早々に吉田知子の本を出版しようとして、作家に「あなたは財閥の御曹司かなにかか」と訝られたそうだ。ははは。

Bookcover 水の出会う場所 [a]
魚住陽子 / 駒草出版 / 2014-09-26

フィクション - 読了:「ポンス・ピラトほか」「水の出会う場所」「わたしの名は赤」

Bookcover 刑務所しか居場所がない人たち : 学校では教えてくれない、障害と犯罪の話 [a]
山本 譲司 / 大月書店 / 2018-05-17
「でも、断言しよう。[知的]障害のある人たちだって、かならずコミュニケーションができる。[...] いちばん共感しやすいのは「悲しい」という感情だ。障害のある人は、悲しい思いをたくさん経験しているから、他人の悲しみに敏感だ。[...] とくに、お母さんが悲しんでいる気持ちはすぐに察する。お母さんの表情が少しでも曇れば、彼らも落ち着かなくなる。そんなとき「悲しいね」「つらいね」って声をかけていると、いつしか心を開いてくれる」

Bookcover コラボ゠対独協力者の粛清 (文庫クセジュ) [a]
マルク・ベルジェール / 白水社 / 2019-11-09

Bookcover 「歴史認識」とは何か - 対立の構図を超えて (中公新書 2332) [a]
大沼 保昭 / 中央公論新社 / 2015-07-24
話し手の先生(故人)は、村山内閣のときにできたアジア女性基金の理事を務めた人。私もまたあの頃は、国による補償ではなく募金だなんて、と反発していたのだが、このような世の中になってしまったいまにして思えば、なにもしないよりは全然ましだったんだろうなあ、と...

Bookcover マトリ 厚労省麻薬取締官 (新潮新書) [a]
瀬戸 晴海 / 新潮社 / 2020-01-16

Bookcover マタイ受難曲 (ちくま学芸文庫) [a]
礒山 雅 / 筑摩書房 / 2019-12-10
夜中にマタイ受難曲を流しながらちびちびと読んだ。楽曲解説にはいると、音楽学の専門用語が頻出しちんぷんかんぷんなんだけど、それでもそれなりに楽しく読了。

ノンフィクション(2018-) - 読了:「マタイ受難曲」「刑務所しか居場所のない人たち」「コラボ=対独協力者の粛正」「『歴史認識』とはなにか」「マトリ」

2020年2月16日 (日)

Bookcover 六四と一九八九:習近平帝国とどう向き合うのか [a]
/ 白水社 / 2019-12-25
昨年に明治大で開かれたシンポジウムを基にした論文集だそうだ。張博樹という民主派知識人(現在は米在住) の章が興味深かった。

Bookcover ハンナ・アーレント (ちくま新書) [a]
森分 大輔 / 筑摩書房 / 2019-06-06

Bookcover ジャパン・ストーリー 昭和・平成の日本政治見聞録 [a]
ジェラルド・L・カーティス / 日経BP / 2019-05-23

Bookcover 新編 戦後翻訳風雲録 (大人の本棚) [a]
宮田 昇 / みすず書房 / 2007-06-02
著者は早川書房の編集者出身。田村隆一がいかに迷惑な人であったかという話が縷々述べられていて、こういっちゃ申し訳ないけど、とても楽しく読んだ。
 鮎川信夫が亡くなったあと、独身だったはずの鮎川の妻として年配の女性(最所フミ)が名乗りでて、周囲は仰天したのだそうだ。それって、ちくま文庫「英語類義語活用辞典」の最所フミではありませんか...

Bookcover 貞観政要 (ちくま学芸文庫) [a]
呉 兢 / 筑摩書房 / 2015-09-09

ノンフィクション(2018-) - 読了:「六四と一九八九 習近平帝国とどう向き合うのか」「ハンナ・アーレント」「ジャパン・ストーリー」「新編 戦後翻訳風雲録」「貞観政要」

<< 読了:「戦争は女の顔をしていない」「ハコヅメ」「あせとせっけん」「きちじつごよみ」「あさひなぐ」
 
validate this page / CSS