{glba}パッケージ備忘録 - 論理の流刑地

体調管理が、何よりいちばん難しいタスク。

Introduction
LBAの概要かいつまみ
Main Function

Introduction

response time（回答時間）をどうやって行動予測や（関心のある変数の）効果の推定精度につなげるか、というハナシをやっている界隈がある。
ご多分に漏れず、この分野にもいくつかの流派があるのだが、そのうちLBA(線形弾道蓄積)モデルを使った分析を提供しているのが{glba}パッケージである。

このパッケージにはよくあるvignette（解説文）がないし、そもそも日本語での情報がほぼほぼないので、色々使いながら試していく。

◆参考URL

RStanによる反応時間の解析：Linear Ballistic Accumulator modelを用いて(専修大：国里先生)

LBAの概要かいつまみ

（本項の内容は、参考URLの国里先生のほうがより端的に、より分かりやすくまとまっているのでそちらも参照されたい。まぁ自分の言葉で書いたほうが定着しやすい、というのをより重視して書いているだけである）

上述のBrown and Heathcote(2008)を一読して大まかな概要を理解した。
簡単に言ってしまえば、複数の回答選択肢*1に対する選択式問題の回答がある一定時間でなされるメカニズムを、「証拠蓄積evidence accumulate」という理論モデルで捉えることにその要諦がある。

問題文や選択課題が提示されたときから回答を行うまでの時間は大別して、

問題・課題の意味を理解するためのインプット刺激の処理のための時間（非決定時間 non decision time,以下non-DTと略記）
問題の意味を解してからの、判断・選択までの時間（決定時間decision time、以下DTと略記）

の二つに分けられる。

LBA（「線形弾道蓄積」）というモデルの名称が意味をもってくるのは、DTの決定メカニズムの表現においてである。
回答時間と回答選択の同時モデリングという課題においては、「どの選択肢が回答者によってえらばれるか」だけが解くべき単一の問いとはならず、同時に「同じ選択肢を選んだ回答者のなかでも回答時間に差異が出るのはなぜか」という問いに対しても答えを用意せねばならない。

そこで、「ある選択肢をある時間で選ぶ」という事象を「非決定時間が終了し決定時間に各回答者は証拠蓄積のモデルで捉え、その個人差を

蓄積のスピード（＝直線の傾き）
蓄積開始地点の値（＝直線の切片）
閾値（＝回答基準）

という三つのパラメータで把捉するのが、LBAだ。
また、各パラメータを変数ではなく（共変量と係数の線形結合としての）関数として推定できるのも、LBAのアピールポイントである

色々すっ飛ばして*2、CDF（累積密度関数）を導出すると、以下のような形になる。

（上記URL1：国里先生の記事より引用）

ここで、式中のnotationは以下の通りである

Aは、直線の切片aのドロー元である一様分布の上限値
bは、閾値
vは、直線の傾きdのドロー元である正規分布の平均値（sは正規分布の誤差の標準偏差だけど、scalingの問題で1に固定

国里先生の至極丁寧な記事では、本モデルでのStanでのbayes実装が公開されていて、そちらも尤度（密度）関数をスクラッチで組んでいてとてもとても参考になる。
しかしここでは、最尤推定（正確にはFIML：完全情報最尤推定法）を用いた推定パッケージとしてRで公開されている{glba}を使ってみる

Main Function

{glba}パッケージに含まれている関数は、

lba()
rlba()
tablba()

の三つで、rlbaは架空データの生成用関数、tablba()はlba()の推定結果を受け取って表示するだけの関数なので、実質{glba}はlba()という一つの関数を提供しているパッケージである。

ひとまず、何はともあれコード例をみてみよう

# loading data
data( "numpp1")
dfNpp <-  numpp1

# setting  starting values and estimation
pars <- c(.25,.3,.1,.2,.9,.75,.6)
m1 <- lba(rt=rt,response=acc,drift=~diff-1,data=dfBH,
          startpars=pars, hessian =T)
#show estimates
tablba(m1)

引数は、以下のように分けられる

データ指定系
- rt：回答時間変数の変数名(character)
- response：回答値変数の変数名(character)
- data ：分析対象となるデータフレーム(data.frame)
- weights ：case weights
モデル指定系（指定しないと単なる関数ではなく単なるパラメータとして各値を推定）
- sddr：drift rateの標準偏差の決定モデル式(formula)
- sp：弾道蓄積の切片の決定モデル式(formula)
- bound：閾値の決定モデル式(formula)
- nond：被決定時間の決定モデル式(formula)
- drift：drift rateの平均値の決定モデル式(formula)
- scaling：識別不能を回避するためのパラメータ制約（"sum" or "fixedSD"）。"fixedSD"を指定した場合、startparsのはじめに標準偏差の値を指定する
- loglink：各モデル（sddr, sp , nond ,bound）にlog link functionを用いるかどうか(logical)
- sdstart：drift rateの初期値。scaling = "fixedSD"のときに使用（startparsの最初の値と同じ値を指定するということ？だと思う）
推定手続き指定系
- startpars：optim()に内部で引き渡すパラメータ推定の初期値(numeric)
- fixed：もしstartparsに引き渡した初期値のまま固定させたいパラメータがあるときはここで指定(logical)
- method：optimの引数と一緒。どの最適化法を使うかの指定（"Nelder-Mead"とか"BFGS"とか）
- hessian：ヘシアンをoutputするか否か(logical)

lba()は内部でoptimを使って最適化していて、そのためパラメータの初期値をstartparsでそのまま指定しなければならない*3。

だから、とくに各パラメータに関して共変量をかませた式を推定するときは、パラメータの数を正しく把捉しなければならない。
このパッケージをたくさん使う必要があるのなら、自動的にstartparsを生成するwrapper functionを作った方がいいかもしれない。
（あと、factor型の変数の扱いがイマイチわかりづらい*4ので、カテゴリ変数はあらかじめ{0,1}のダミー変数化して投入すべきかもしれない）

上記コードの、出力結果は以下のようになる