コムアーチのショールカラーカーディガン。
横浜仲町台のEuphonicaさまで購入。
東北地方にて昔ながらの手横編み機のみでニット製造を行うニッターさんの手で丁寧に編み立てられて作られたものだとか。
襟部分の返りが柔らかく立体的で、首元を優しく包んでくれる。羽織るとずっしり重さを感じるものの、不思議と負担感はない。
あとカーディガンとは思えない威圧感とオーラがある。
一流のニットデザイナーと一流の職人が作り上げた逸品。
良いものには全て根拠があるのだ。
コムアーチのショールカラーカーディガン。
横浜仲町台のEuphonicaさまで購入。
東北地方にて昔ながらの手横編み機のみでニット製造を行うニッターさんの手で丁寧に編み立てられて作られたものだとか。
襟部分の返りが柔らかく立体的で、首元を優しく包んでくれる。羽織るとずっしり重さを感じるものの、不思議と負担感はない。
あとカーディガンとは思えない威圧感とオーラがある。
一流のニットデザイナーと一流の職人が作り上げた逸品。
良いものには全て根拠があるのだ。
昨今有名な話ではありますが、虐待によって子どもの脳は変形します。
https://toyokeizai.net/articles/-/189445
脳が変形するということは、本来有している脳の機能がうまく発揮できず、本来できるはずのことができなくなるなど、様々な困難が生じてしまうということです。
以下、虐待と脳についての分野で研究を続けられている友田明美先生のコメント(抜粋)です。
虐待やネグレクト(育児放棄)などの不適切な養育は愛着障害を引き起こします。
日常的に養育者が子どもに暴言虐待や長期的な厳格体罰を与える、そうすると不安定な愛着が形成されてしまいます。
~
養育者が戻って来たときも同様で、子どもは喜びもしないどころか、そっぽを向いたままです。このように不適切な養育が引き起こす愛着障害は、こころの発達に問題を抱え、さまざまな症状を表します。
症状が内向きに出ると、他人に対して無関心になったり、用心深くなったり、イライラしやすくなったりして、他人との安定した関係が築けません。症状が外向きに出ると、多動で落ち着きがなくなったり、友達とのトラブルが多くけんかが絶えません。また、礼儀知らずとなり、対人関係に支障をきたしてしまいます。実は、虐待などが原因で、社会的養護を受けている子どもたちの40パーセントに愛着障害が発症することがわかってきました。
幼児期に虐待ストレスを受け続けると、脳の中にある感情の中枢である扁桃体(へんとうたい)が異常に興奮し、副腎皮質にストレスホルモンを出すよう指令を出すのです。
そうするとストレスホルモンが過剰に放出され、脳にダメージを与えるのです。
アメリカのハーバード大学との共同研究でわかってきたことは、感情をつかさどる前頭葉が小さくなって、自分のコントロールができなくなり、凶暴になったり、集中力が低下したりします。暴言虐待により聴覚野が変形し、聴こえや会話、コミュニケーションがうまくできなくなったりします。両親間のDV・家庭内暴力を目撃すると視覚野が小さくなり、他人の表情が分かりにくくなり、対人関係がうまくいかなくなったりします。脳の形が変わるのは、「外部からのストレスに耐えられるように情報量を減らす」ための脳の防衛反応だと考えられています。
http://www.nhk.or.jp/kaisetsu-blog/400/264452.html
このように、様々な脳部位に影響を及ぼすのが虐待です。友田先生の研究では他に、脳の感受性期の関係で、被虐待の時期によりダメージを負う脳部位が異なることも示されています。
Executive function performance and trauma exposure in a community sample of children, Anne P.DePrincea.etでは、虐待と、ワーキングメモリ、抑制、聴覚注意、および処理速度タスクで構成される実行機能のパフォーマンス低下との関連が示されています。
https://www.sciencedirect.com/science/article/abs/pii/S0145213409000969
Cognitive impairment in school-aged children with early trauma,JoanaBücker.etでは、wiscの数唱の短さ(ワーキングメモリの下位検査)が指摘されています。
https://www.sciencedirect.com/science/article/abs/pii/S0010440X11002392
一方で、こちらでは言語理解と処理速度が比較群より低く、知覚推理とWMは同レベルという結果が出ています(COGNITIVE ABILITIES OF MALTREATED CHILDREN
,Kathleen D. Viezel Benjamin D. Freer Ari Lowell Jenean A. Castillo)。
https://onlinelibrary.wiley.com/doi/abs/10.1002/pits.21809
いずれにせよ、虐待により認知機能の低下は避けられず、虐待の内容や被虐待年齢等により脳の器質的変化が多様であるという友田先生の有名な研究より、認知機能の変化も多様ということになるのかもしれません。
今更なんでこんな話題をっていうと、
とあるツイートで興味深いものを拝見したからでした。
https://twitter.com/tsuyomiyakawa/status/1217400388724248576
脳トレ的な認知機能トレーニングに関する多数の研究&メタ解析の結果、その種のものには一般的認知機能を向上させる効果はほぼ認められないというものでした。
その後のツイートの中で、発達障害も含め、各種の脳の疾患においての一般的な認知機能に対する認知機能トレーニングの効果は、この総説では検証されていないと続けられています。
虐待により脳に直接与えられたダメージを回復するには、こういった認知機能を回復させる認知機能トレーニングが有効じゃないかと感じており、高次脳機能障害者への認知機能トレによるエビデンスを参考にしていたところでしたので、
一瞬この研究にはびびりましたが、後半のコメントによりまだ可能性が残されていると思い少し安心。
現状から機能向上は難しくても、本来有している機能が損なわれている状態から元の状態に回復する。そういうイメージであれば、認知機能トレーニングが有効であればいいなと思いますし、
偉い人がここらへんを検証して、児童虐待によるダメージから回復するエビデンスが構築されていくと児童虐待の世界は飛躍的に進歩するのになぁ、と思います。
ここらの認知機能が改善したら、情動ラベリングや行動の心理的意味の言語化とか、前頭葉機能を用いる被虐待により歪んだ認知改善もスムーズになるんですが、それはまた別のお話。
項目反応理論(IRT)による分析のモデルの1つに段階反応モデルがあります。
「はい」「いいえ」のような2値反応ではなく、「はい」「どちらでもない」「いいえ」のような3つ以上の順序尺度的な反応を求められる課題があるとします。それをIRTの枠組みの中で分析する方法が、段階反応モデルになります。
やや小難しい説明になりますが
段階反応モデルでは、項目 j(=1, …, J ; J は項目数)の反応ujがC個の値をとる順序尺度の離散変数であると仮定します。
以下に、課題a~eの得点を、3つの順序尺度で得点化したデータがあるとします。たとえば、「a. 算数より国語の方が好きだ いいえ=0、どちらでもない=1、はい=2」「b. 算数より国語の方が高得点の傾向にある いいえ=0、どちらでもない=1、はい=2」のような質問紙があるとして、以下のように記入してもらうという感じです。
a b c d e
A君 1 2 3 2 3
B君 3 3 3 2 3
C君 1 2 1 1 2
D君 1 1 2 2 1
E君 3 1 3 2 3
F君 2 1 3 2 3
G君 3 3 3 2 3
H君 1 3 3 2 3
・
・
・
このデータを、段階反応モデルを用いた項目反応理論による分析を行うとします。
結果と解釈は以下のようになります。
各課題を、いいえ=0、どちらでもない=1、はい=2として3カテゴリの順序尺度とし、MCMCによるIRT段階反応モデルのパラメータ推定を実施。
、
その結果
課題4は、困難度の観点で見ると、b「4,1」とb「4,2」の差が最も大きい。これは中間に位置するカテゴリ(1点)に反応する確率が最も高いことを示しており、すなわち偶然正答の「1」の可能性の高さが示されている困難度の差が最も大きく「どちらでもない」の可能性が最も高い課題であることが示された。
課題1は適切な困難度かつ、偶然正答の「1」の可能性も最も低く、個人特性を最も識別する課題であることが示された。
課題3、課題5については、困難度の最下位(b[3,1]b[5,1])、最上位(b[3,2]b[5,2])が負の域に入っていることから、困難度の低い課題、すなわち平均以下の個人特性でも正答する率が高い傾向にある課題であることが示された。
各課題の特性を検討するには項目反応理論は有用で、段階反応モデルはその拡張版といえます。
個人的には、検査をある種質的な観点で解釈できる分析手法として好んで用いています。でも使いこなせているかはまた別の話…。
本解析で用いたStanコードはこちら。
data{
int ni;
int nj;
int nc;
real D;
int<lower=1,upper=3> y[ni,nj];
}
parameters{
vector[nj] a;//識別力母数
ordered[nc-1] ba[nj];//困難度母数
vector[ni] theta;//回答者の特性値
}
transformed parameters{//項目反応カテゴリ特性曲線の定義を用いる
real b[nj,nc];
vector<lower=0,upper=1>[nc-1] pa[ni,nj];
simplex[nc] p[ni,nj];
for (j in 1:nj){
for (c in 1:nc){
if (c ==1){//もしcが1(最下位カテゴリなら)
b[j,c]<-ba[j,c];
}else if (c ==nc){//cがnc(最上位カテゴリ)なら
b[j,c]<-ba[j,c-1];
}else{//cが上2つ以外なら
b[j,c]<-(ba[j,c-1]+ba[j,c])/2;
}
}
}
for (i in 1:ni){//i:1~ni
for (j in 1:nj){//j:1~nj
for (c in 1:nc-1){//c:1~nc-1
pa[i,j,c]<- 1/(1+exp(-D*a[j]*(theta[i]-ba[j,c])));
}
}
}
for (i in 1:ni){
for (j in 1:nj){
for(c in 1:nc){
if (c==1){
p[i,j,c]<-1-pa[i,j,c];
}else if(c==nc){
p[i,j,c]<- pa[i,j,c-1];
}else{
p[i,j,c]<- pa[i,j,c-1]-pa[i,j,c];
}
}
}
}
}
model{
for (i in 1:ni){
theta[i]~normal(0,1);
for (j in 1:nj){
y[i,j]~categorical(p[i,j]);//カテゴリカル分布:結果が3種類以上の状態を持つ施行が従う分布
}
}
for (j in 1:nj){
a[j]~lognormal(0,sqrt(0.5));
for (c in 1:nc-1){
ba[j,c]~normal(0,2);
}
このとき、一般的な項目反応理論の説明と同様に、潜在能力θの被験者がuj =cと反応する確率Pjc (θ)は、潜在能力θの被験者が uj≧cと反応する確率Pjc(θ) を用いて表現できます。
θを変数とみると段階反応モデルにおける項目jのカテゴリcの項目特性曲線を表しており、これは項目反応カテゴリ特性曲線(Item Response Category Characteristic Curve; IRCCC)と呼ばれます。
また、境界特性曲線(Boundary Characteristic Curve; BCC)は、θの値によらず
Pj0(θ)=1 ・・・0より大きくなるのは確率1(100%)
PjC(θ)=0 ・・・Cより大きくなるのは確率0(0%)
です。
ポールスミスと聞いて浮かぶイメージとは、やはり派手な裏地のジャケットやカジュアルに振ったマルチストライプの財布でしょうか。
これはコレクションラインの、いわゆる総柄プリントシャツ。
コレクションラインの良さは、多くはメイド・イン・ジャパンで、素材の良さや縫製・着心地の良さなどもよく語られますが
個人的には、Paulsmith Collectionの持つ退廃的な世界観だろうと思います。
構造方程式モデリング(SEM)は、多変量解析と呼ばれる統計手法群の1つですが、この手法の特徴は柔軟なモデル構成力にあります。
ある概念を測定したいと心理検査を構成した際に、その概念の下位にある要素間の関係性や概念に対して要素がどのように影響しているか、などを確認する手法です。
SEMの特徴としては他に、「構成概念f」と「観測変数x」があります。観測変数というのは実際の尺度項目、構成概念は実際にある項目ではないですが観測変数により確認された仮想的な変数です。
数学、英語、国語という科目が観測変数、学力が構成概念、という感じです。
SEMは研究論文などでなんとなく見たことある方も多いとは思いますが、
といったところは正直よく分からず見ている方が多いのではないでしょうか。
せっかくなので、以下の図を例に解説していければと思います。
図(考え方を簡単にするため簡略化)
e:独自変数…観測変数の固有にもっている変数で、構成概念ではカバーできない範囲
v:観測変数…実際に測定したもの
f:構成概念…観測変数に影響を与えている概念
a:影響指標・因子負荷・パス係数etc…構成概念から観測変数への影響力の強さ
この図の場合、3つの測定方程式が立てられます。
v1=a1*f1+e1
v2=a2*f1+e2
v3=a3*f1+e3
と表現されます。モデルを分かりやすくするため、(本来はそんなことないんですけど)基本的な「学力」で説明しきれない個々のテストの成績は互いに無相関、ということにして、
E[eiej]=0 (i≠j)
と仮定します。
観測変数の共分散と測定方程式の母数との関係を記述すると(途中式省略)
σ12=E[v12]=a12+σe12
と表現されます。同様に
σ22=a22+σe22
σ32=a32+σe32
となります。次にv1とv2の共分散は(途中式省略)
σ21=E[v2v1]=a2a1
σ31=a3a1
σ32=a3a2
が導かれます。
観測変数の共分散行列の要素がこれで全て揃いました。以下に整理すると
と構造化して書き直すことができます。上のように共分散を方程式モデルの母数で表現したものを「共分散構造」といいます。
この共分散構造行列がSEMの係数導出等の基本になるものでして、具体的な集計データからこの共分散構造行列にあてはめ、観測変数の標本共分散行列を構成して、最終的な係数等があてはめられたモデルが作られていくという感じです。
次に、パス係数や因子負荷などの推定量の導出です。
a1=±√(σ21σ31)/σ32 が得られます。a1が導出できれば残りの係数も導出できます。
a2=σ21/a1
a3=σ31/a1
誤差分散は
σe12=σ12-a12
σe22=σ22-a22
σe32=σ32-a32
より導出されます。
その後、分散が1に標準化された「標準化モデル」により、モデルを解釈しやすい形に数値を変換する、などの作業を行ったりします。
以上のようにSEMは、潜在した共通原因のために観測変数の間に相関が生じる現象を表現した「測定方程式モデル=因子分析モデル」から共分散行列を構成し、その共分散行列を利用して測定方程式モデルの母数の推定量を導出することで、各観測変数や潜在変数(構成概念)の関係性をモデルで表現するものといえます。
SEMのパスとかって一体どこから出てくるんだろうとか思っていた時期がながくありましたが、この共分散行列に由来していたんですね。
こうやって作られたモデルの適合度を測定する指標は色々ありますが
SEMで使用するものとしてはRMSEAとか、我らが赤池先生の考案したAICとかが有名でしょうか。
RMSEAはSEMに特化して、モデルの分布と真の分布との乖離を、1自由度当たりの量として表現した指標になります。
RMSEA=√max((fML/df)-(1/N-1),0) ※fML=最尤法の目的関数、df=自由度
0.05以下であればあてはまりが良く、0.1以上であれば当てはまりが悪いとすることが多いように思われますが、この数値自体に根拠があるわけではないので注意。
AIC=χ2-2df
尤度で定義された統計モデルの良さを測定するために使用され、数値が小さいほど良いモデルだと判断されます
なお、χ2=(N-1)fMLで定義され、適合が悪いほど値が大きくなります。カイ自乗値(検定統計量)はより小さく、p値はより大きなものであることが望ましいといえます。値そのものを適合度の指標として用いるよりは、「モデルはデータに適合している」という帰無仮説を検定するための検定統計量として用いられます。サンプルサイズが大きければ得られるp値も大きくなります(つまり第1種の誤りの確率が高くなる)が、SEMではそもそもサンプルサイズの大きいデータを用いることが多いので、この検定結果に意味を求める感じではないです。
SEMの適合度指標については、星野・岡田・前田(2005)『構造方程式モデリングの適合度指標とモデル改善について:展望とシミュレーション研究による新たな知見』が詳しいのでそちらを参照に!
因子分析には様々な技法があります。
1個人に連続してデータをとって時系列データっぽいものを対象としたp技法因子分析とか、サンプルが集まらない際の苦肉の策として私個人としても用いたことがある技法もあったりしますが、
今回は「延長因子分析(extention factor analysis)」のお話をしたいと思います。
大体が『延長因子分析の方法論― 変数と因子との相関係数として定義される因子構造を用いて―』清水和秋(2012),関西大学心理学研究,2012,第3号 からの情報なので、ここを見た方が正確ではあるんですが、せっかく論文を読んで学んだことのまとめということで。
そもそも延長因子分析って、実務的には追加した変数をパス関係に置くことによる測定モデルの意味とか、イマイチ腑に落ちない部分もある(理解不足)ので、延長因子分析を用いるメリットデメリットなどをどうぞアホな私に教えて偉い人。
我が恩師の清水先生の論文は無駄がなく簡潔なので、自分のようなアホにはようわからんことが多いのです。
延長因子分析って、元々行われた因子分析結果があったとして、そこに新しい変数を加え、その新しく加えた変数の因子パターンをどうやって分析するか~っていうものなんです。
その分析の方法としていくつか方法論が提案されているんですが、
2000年前後から広まってきた構造方程式モデリング(=SEM)を用いるとかなり簡便に延長因子分析が可能となる、
そんな感じでよろしいのでしょうか偉い人。
因子分析の基本モデルを行列で表すところから始めます。
まず、あるN人の被験者について、n変数の測定を行い、m 個の因子を抽出することができたものとします。次に、この観測変数の標準得点行列をZ(N×n)とすると因子分析のモデルは次のように表すことができます。
F は(N×m)次の因子得点行列、
Vfpは(n×m)次の因子パターン行列、
(1)式のaj1~ajm を第j 番目の変数の因子パターンm 個の値とすると、
この行列の要素(n×m)はn 個の変数の因子パターンからなります。
この共通因子空間とは独立した独自性に関するものが(N×n)次の独自性得点行列U と、独自性を対角項にもつ(n×n)次の対角行列Dになります。
R を(n×n)次の変数間の相関行列とすると、次のように展開することができます。
= VfpCfVfp +D2
Cf は(m×m)次の因子間相関行列です。
(2)式は、因子得点間の相関行列でもあります。
Vfp は初期の因子解を回転した結果の因子パターン行列です。
これに対して、変数と因子との相関係数が因子構造であり、この行列Vfs(n×m)(=因子構造行列)は、次のように計算することができます。
・(3) Vfs =(1/N) Z'F
因子パターン行列と因子構造行列との関係は次のように表すことができます。
・(4)Vfs = Vfp Cf ⇒ 因子構造行列=因子パターン行列×因子間行列
ここでは、この得点の推定値Fˆ を算出したものとします。そして、同じN 人を対象としてi 個の変数を追加することができたものとします。
この標準得点行列をiZ(N×i)と表して、上の因子得点との相関係数を(3)式のように求めてみると以下のようになります。
・(5)iVfs =(1/N)Z Fˆ
この(5)式で得られた新しい変数の因子構造行列に(4)式のように因子間相関行列の逆行列を次のように掛けると
・(6) iVfp= iVfs C-1f
となりまして、新しく追加した変数の因子パターンの値iVfpを得ることができました。
すなわち、因子分析での対象ではなかった新たに加えられたi 個の変数をm 次元のn 変数から抽出された共通因子空間に布置させることができた、ということになります。
以上が、古典的な延長因子分析の考え方になります。
◆SEMによる方法 ※図は清水(2012)から引用
イメージとしては、以下の図のような感じです。追加変数を設定して、まとめてSEMによる分析を行います。
Fig.1 8 変数2 因子の測定モデルに3 つの変数を追加
標準化解:因子パターンp1 ~ p8、因子間相関
延長因子分析:因子と追加した変数との相関係数r1 ~ r3
図のr1 ~ r3 の因子と追加した変数との相関係数は、(5)式のiVfs(iの因子構造行列) に相当します。これらの値は、(6)式のように因子パターンへの変換が必要となります。
※直接iVfp (新しく追加した因子パターンの値)を推定するために測定モデルである因子からのパスとしてこれらの追加した変数を置くことも可能なように思えますが、このやり方では因子を構成する変数に変化を起こすことになるそうなので、ここでは追加した変数は相関関係としておいています。
SEM の方法からは、このように、因子得点と追加した変数との相関係数を直接的に計算することができ、この点が古典的な方法とは異なる点になります。適合度等の検定も可能となるので、延長因子分析についてはSEMを用いると便利なように思えます。
昔の偉い人は色々工夫をして繁雑な計算を行っていたのですが、複雑な計算を行わずに済むようなソフトを開発してくれた最近の偉い人のお陰で自分のようなアホに様々な恩恵があると思うと、今と昔の偉い人はすごいなと思います。
児童虐待対応の仕事、いわゆるケースワーカー(CW)の仕事は日々とても辛いものがあります。
でもどうして辛いんでしょうか。
仕事の内容の激しさももちろんありますが、その他のとてもシンプルな理由に、仕事の多さ、次々と終わらない担当の受理などがあります。
要するに、この前の仕事が終わっていない(すぐ終わるものでもないし)のに次々くる感じです。
これはもう誰でも嫌になっちゃいますよね。
そんなシンプルなつらさってどうやってモデルで表現するのでしょうか。
今回、色々と小難しいことは抜きにした「CW嫌々モデル」を考えてみました。
架空データとして、「嫌だなって気持ちを点数化したもの(高いほど嫌)」である嫌得点(iyapoint)と、仕事が増えた・殺害予告された・個人名で訴えると怒鳴られた・殴られた(どれも日々のあるあるだね!!)のイベントが発生したかを示す2値データ(ivent)が、日々の時系列データでどのように変化するかを見ていきます。
day iyapoint ivent
1 2019/4/1 7 0
2 2019/4/2 7 0
3 2019/4/3 9 1
4 2019/4/4 8 0
5 2019/4/5 9 1
6 2019/4/6 9 0
7 2019/4/7 8 0
↑こんな感じのデータです。
早速モデルを考えていきます。
◆モデル1:トレンド+嫌イベントストレスモデル
1.放っておくと嫌々度は下がる;-Trend(負の傾きの直線をイメージしてください)2.定期的に負のイベントがある:CoefIvent*Iventt (トレンドの係数=一定)×嫌なイベントの有無
1,2よりイベントの影響を含めたトレンド:δt=Trend+CoefIvent*Iventt
3.μt=μt-1+δt-1+wt, wt~(0,σ2w):μtはt時点の嫌々度の状態
4.μnoiset=μt+Vt:Yt,Vtは各々t時点の過程誤差と観測誤差
5.λt=exp(μt)
6.Yt~Poisson(λt)
嫌なイベントの発生はまれ(だといいなという想定)という前提で、ポアソン分布を使用。でも嫌々は得点化したから、ポアソンのような離散分布でなく連続データの分布にした方がよかったかなぁと今更ながら思った。
少し雑念が入りましたが、上記1~6を元に最初の数理モデルを組み、統計処理をしてみようと思います。
統計ソフト:R-3.5.2
使用パッケージ:Rstan他
連続データを用いた、状態空間モデルによるベイズ推定を行いました。
以下Stanコードを記します。
stancode <- "
data {
int n; // サンプルサイズ
int ivent[n]; // 嫌イベントの有無
int iyapoint[n]; // 嫌得点
}
parameters {
real trend; // 嫌得点減少トレンド
real coef_ivent; // 嫌イベント発生がトレンドに与える影響
real mu_zero; // 状態の初期値
real mu[n]; // 状態の推定値
real mu_noise[n]; // 観測誤差の入った状態の推定値
real<lower=0> s_w; // 過程誤差の分散
real<lower=0> s_v; // 観測誤差の分散
}
transformed parameters{
real delta[n]; // イベントの影響が入ったトレンド
real lambda[n]; // ポアソン分布の期待値
for(i in 1:n){
delta[i] = -trend + coef_ivent * ivent[i];
}
for(i in 1:n){
lambda[i] = exp(mu_noise[i]);
}
}
model {
// 状態の初期値から最初の時点の状態が得られる
mu[1] ~ normal(mu_zero, sqrt(s_w));
// 状態方程式に従い、状態が遷移する
for (i in 2:n){
mu[i] ~ normal(mu[i - 1] + delta[i - 1], sqrt(s_w));
}
// 観測誤差が加わる
for(i in 1:n){
mu_noise[i] ~ normal(mu[i], sqrt(s_v));
}
// ポアソン分布に従って観測値が得られる
for(i in 1:n){
iyapoint[i] ~ poisson(lambda[i]);
}
}
generated quantities{ //新たにサンプリングする変数を作る
vector[n] log_lik;
for(i in 1:n){
//モデル1の対数尤度
log_lik[i] = poisson_lpmf(iyapoint[i] | lambda[i]);//ポアソンの場合、lpdfではなくてlpmfであることに注意。https://www.slideshare.net/simizu706/stan-64926504の42
}
}
"
fit <- stan(model_code = stancode, iter = 8000, chains = 1, data = standata, seed=123)
◆モデル2:CW嫌イベント確率的ボラティリティ変動モデル
確率的ボラティリティモデルとは、金融界でよく使われるモデルです。ボラティリティ(分散的な)が常に一定であるランダムウォークとは異なり、ボラティリティ自体もランダムウォークするというモデルです。株価みたいに、上がる時はじわじわ上昇トレンドで上がるくせに、下がる時は理不尽にガタッとくる値動きをよく表現できているモデルだと思います。
何もなければ(ありえないけど)すこーしずつストレスが和らぐけど、ちょっとしたことで一気に限界突破する児童虐待の仕事と似てると思いませんか?
以下、モデル2のStanコード
stancode2 <- "
data {
int n;
real iyapoint[n];
}
parameters {
real h[n];
real mu;
real<lower=-0.999,upper=0.999> phi;//ボラティリティの持続性
real<lower=0.0001> sigma;
}
model {
h[1] ~ normal(mu,sigma/sqrt(1-phi*phi));
iyapoint[1] ~ normal(0,exp(h[1]/2));
for (i in 2:n){
h[i] ~ normal*1;
}
}
generated quantities{ //新たにサンプリングする変数を作る
vector[n] log_lik;
for(i in 1:n){
//モデル2の対数尤度
log_lik[i] = normal_lpdf(iyapoint[i] | 0,exp(h[i]/2));
}
}
"
fit2 <- stan(model_code = stancode2, iter = 8000, chains = 1, data = standata, seed=1234)
そして、このモデル1とモデル2のどちらがデータに適合したモデルなのでしょうか。
適したモデルこそ、児童虐待CWの辛さをよく表現したモデルってことになります。
モデルの比較はWAICを使用します。
適合度比較では、構造方程式モデリングではAICとかGFIとかRMSEAとかがスタンダードですが、ベイズではWAICがよく使用されています。以下にRで実行する際のコードを示します。
log_lik1 <- extract_log_lik(fit)#waic用コードで計算したfit12の対数尤度の値を取り出す
waic1 <- waic(log_lik1)#変数名を付ける
print(waic1 , digits = 3)#waic1の結果表示
log_lik2 <- extract_log_lik(fit2)
waic2 <- waic(log_lik2)
print(waic2 , digits = 3)
compare(waic1,waic2)#waicの比較.昇順に並べ替える
print(waic1 , digits = 3)#waic1の結果表示
Estimate SE
elpd_waic -198.555 1.014
p_waic 1.131 0.087
waic 397.111 2.028
log_lik2 <- extract_log_lik(fit2)
waic2 <- waic(log_lik2)
print(waic2 , digits = 3) #waic2の結果表示
Estimate SE
elpd_waic -337.736 1.597
p_waic 0.161 0.013
waic 675.471 3.195
compare(waic1,waic2)#waicの比較
elpd_diff se
-139.2 0.9
★結論:モデル1(トレンド+嫌イベントストレスモデル)のWAICの方が小さいため、予測の観点から良いモデルといえることが分かりました。
とはいえ、モデル2の確率的ボラティリティな感じも、日常業務である「落ち着いてきていけると思ったけど炎上してやっぱ無理」な感じをよく表現できていると思います。
確率的ボラティリティを用いた、より良いモデルに改良していけるのではないかと思います。今回のはシンプルでしたしね。