筋トレと状態空間モデルと私



この記事は Stan Advent Calendar 2017 の12月22日の記事です。私の筋トレのデータをつかって状態空間モデルで遊んでみた、という内容です。今年は状態空間モデルの素晴らしい記事が多いので恐縮ですが。すみません、僕には筋肉しかないんです。

基本的にHijiyamaR#finalでLTしたものと同じです。#ありがとうぞうさん

 

筋トレと私

体を動かすことが好きで、院生の時は走ったりしていましたが、あまり時間が取れなくなり、短時間でちゃんと運動した感が得られる筋トレをしていました。最近はあまりトレーニングできていません。好きな種目はデッドリフトです。

 

状態空間モデルと私

広島ベイズ塾の2016年の夏合宿で、山口大学の小杉先生から状態空間モデルを教えてもらいました。昨年のStan Advent Calendarに小杉先生の書かれたStanで体重の推移をみつめてみた(状態空間モデル)という記事があります。面白いなと思って、それからトレーニングデータを記録しておくことにしました。今回あつかうデータは、デッドリフトの第一セットの使用重量と挙上回数の67時点のデータです。使用重量と挙上回数から最大挙上重量(1RM)が計算できるので、そのデータを使います。デッドリフトをするのはよくて週2回とかなのでデータがたまるまで結構かかってしまいました。あと、日付を記録しておくのを忘れていたため、時点の間隔が一定ではないのが問題です。次から気をつけます。

 

筋トレと状態空間モデル


状態空間モデルは今年のStan Advent Calendarで何回もでてきていますので、詳しくはそれらの記事を参照してください。
状態と観測値をわけて、1) ある時点での観測値はその時点での状態に誤差がのってあらわれてくる、2) ある時点の状態は1つ前の時点の状態から変動する、と考えます。


状態空間モデルは筋トレを例にすると説明しやすい。と私は思います。状態=筋力、観測値=パフォーマンス、と考えればよいのです。
最大挙上重量をプロットしたものが下の図です。けっこうガタガタです。筋力がそんなばらついているはずがない。ですよね?


筋トレのパフォーマンスは結構ぶれます。たとえば同じ週でも、前回扱えた重量がめっちゃ重く感じたりします。睡眠不足だったり、疲れてたり、気分がのらなかったり、ご飯食べてなかったり、そんなこんなで、筋力の変動はある程度小さくても、パフォーマンスはわりと変動するのです。知りたいのは筋力の変動なので、パフォーマンスの変動と分けて考えれる状態空間モデルはとても便利。

 

以下では、3つのモデルを扱います。


内容についてはJ.J.F. コマンダーの状態空間時系列分析入門を参考にしております。間違いがあればご指摘ください。


ローカル・レベル・モデル


一番シンプルなモデルです。t = 1,…, nに対して次のようにあらわせます(tは時点、yは観測値、μは状態)。状態のレベル要素が時点によって局所的に変動することを認めるので、ローカル・レベル・モデルといいます。


ローカル線形・トレンド・モデル


ローカル・レベル・モデルに傾き(トレンド)の要素 (νt) を加えたものです。t = 1,…, nに対して次のようにあらわせます。レベルをモデル化するものと傾きをモデル化する式の2つの状態方程式があります。傾きの要素が時点によって局所的に変動することを認めるので、ローカル線形・トレンド・モデルといいます。


グローバル・トレンド・モデル


ローカル線形・トレンド・モデルでは傾きの時点による変動を認めていますが、「傾きは固定」と考えることもできます。t = 1,…, nに対して次のようにあらわせます。時点によって変化しない固定された傾き (ν) を加えています。傾きの要素は時点によって局所的に変動することを認めない(=大域的に決まっていると考える)ので、グローバル・トレンド・モデルと呼んでおきます。なんと呼ぶのが一般的なのでしょうか。

 


これら3つのモデルで分析した結果です。


μtの80%信頼区間と予測分布をグラフにしたものです。グラフの書き方については、アヒル本 StanとRでベイズ統計モデリングを参考にしております。おせわになっております。

左からローカルレベル、ローカルトレンド、グローバルトレンド、で分析したグラフです。予測分布のかたちを見てください。


ローカルレベルでは、前の状態に誤差がのっていくだけなので、予測分布の中央値(黒線)は推定された最後のデータ時点のμの中央値156.58からほとんど動いていません。


これに対して、ローカルトレンドとグローバルトレンドでは、傾きの成分がありますので予測分布の中央値が上がっていっています。ローカルトレンドで推定された、最後のデータ時点のμの中央値は156.28で、νの中央値は0.56でした。グローバルトレンドで推定された、最後のデータ時点のμの中央値は156.38で、(時点で変動することを認めない)固定値としての傾きνの中央値は0.59でした。ローカルレベルでの最後の時点での傾きの推定値とグローバルトレンドでの固定した傾きの推定値がほとんど一緒の値なので、予測分布の中央値は同じような傾きになっています。ただし、ローカルトレンドでは、レベルの変動だけでなく傾きの変動も加わりますので、予測分布の広がりは傾きの変動を仮定しないグローバルトレンドよりもどんどん大きくなっていきます。


ローカルトレンドとグローバルトレンドの違いが見えにくいかと思いますが、これはデッドリフトのデータがたまたま最後の傾きと全体的な傾きが似ていたことによります。ちょうどオーバーヘッドプレスのデータ(87時点)があったので、このデータで分析した結果が下の図です。


オーバーヘッドプレスはデータの最後の方は飽きがきていて、同じ重量で同じ回数こなす作業、という感じでした。そのため、グローバルトレンドではデータの最初の方の伸び率も含まれますので予測分布の中央値は上がっていますが、ローカルトレンドでは最後の方の作業感の影響をうけて傾きがない状態になっています。傾きの推定値の中央値は、グローバルトレンドでは0.21、ローカルトレンドでは-0.03(最後の時点の傾き)となっていました。

 

筋トレと状態空間モデルと私

筋トレのデータを3つの状態空間モデルで分析してみて、モデルの違いが実感できました。ここからは私の個人的な問題になるのですが、筋トレデータに状態空間モデルをあてはめて分析したいと思ったのは、推定された筋力に基づいて、今日のトレーニングの内容を決めたかったからです。


筋肥大に適切な負荷は最大挙上重量の70%から85%(6回から12回ぎりぎり反復できる重さ)といわれています。重さは前回のトレーニングでどれくらいを扱えたのかを参考に「今日はこの重さで何回くらいかな」と決めるのですが、これをもっとちゃんと決めたい。


あとは発生させたMCMCサンプルをつかって計算するだけなので簡単です。ただ、バーベルにつける重りは一番軽くて1.25kgです。なので、2.5kgきざみでしか調整できません。また、細かい重量設定は個人的に結構めんどくさいのです(重りをつけるだけなんですが)。なので、扱いたい重さを決めて、何回やればいいのかを提案してもらうようにしました。たとえばMCMCサンプルの平均値をつかって計算すると、120kgでトレーニングしたい場合には10回を目標にすればよいことがわかります。ただ、せっかく現時点での筋力が分布で得られているので、回数も分布で知りたいですよね。発生させたMCMCサンプルをすべて使って挙上回数を計算し、ヒストグラムにしたのが下の図です。


「普段よりもちょっと負荷をあげる」という漸新性過負荷の原則に従えば、11回を目標にするとよいかもしれません。また、7回以上はほぼ確実に(99.5%)あがることがわかります。仮に、9回しかあがらなかったとしても、その確率は36.4%なので、まぁそんなもんかな、と思うこともできます。

 

まとめ


このように状態空間モデルをベイズでやることによって、今日のトレーニングで何回を目標にするのかが決めやすくなり、また挙上回数の変化に一喜一憂することもなくなるかと思います。ちなみに、今回のような使い方をするのであれば、ローカルレベルモデルで十分かなと思いました。ローカルトレンドモデルだと、特に傾きの変動の標準偏差の収束が悪く、結構時間がかかりました(といっても10分くらいなのですが)。


今回扱ったモデルのコードは続きにおいています。コード等は小杉先生の記事とアヒル本 StanとRでベイズ統計モデリングを参考にさせていただきました。私が手を加えた部分もありますので、間違いがあればご指摘ください。

Continue reading “筋トレと状態空間モデルと私”

「なにかやりたい気持ち」改め、
「説明変数にノイズを含むモデルで遊んでみた」



この記事は Stan Advent Calendar 2016 の12月10日の記事です。
Stan Advent Calendar 2016 盛り上がってますね!!楽しいです!!

アヒル本 StanとRでベイズ統計モデリング の7.7「説明変数にノイズを含む」のモデルで遊んでみた、という内容です。モデルについての誤解がある気がしています。お気づきの点がございましたら、ご指摘いただければ幸いです。
 
 

説明変数にノイズを含むモデル


このモデルは説明変数をX、目的変数をYとしたときに、以下のようになります。
\begin{eqnarray}
X[n] & \sim & Normal(X_{true}[n], 1.0)\\
Y[n] & \sim & Normal(\alpha+\beta X_{true} [n], \sigma)
\end{eqnarray}


はじめの式は「説明変数が平均を真の値(Xtrue) 、標準偏差を1.0とした正規分布から発生している」、つまり「データとして得られた説明変数には標準偏差1.0のノイズが含まれている」ということを表現しています。ここではノイズについて1.0という固定値を与えていますが、もちろん変更できます。
二つ目の式では、説明変数の真の値で目的変数を回帰しています。アルファが切片、ベータが回帰係数、シグマが誤差標準偏差を示します。


構成概念をあつかう心理学では、説明変数に測定誤差などのノイズが含まれることが多いです。そこで複数項目で測定を行い、構造方程式モデリングで誤差と分離して目的変数への影響を検討しているわけなのですが、既存尺度を利用してデータを得た際に因子分析とかを各研究がそれぞれやり直しているような状況です。先行研究の因子分析の結果からは「どの項目にどの程度因子が寄与していて、ノイズがどのくらい含まれるのか」ということがわかるわけですから、そういった情報を利用して「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思っていました。因子分析は潜在変数が説明変数、観測変数が目的変数となっている回帰分析と捉えることができます。はじめの式は「研究者が見たい真の値に誤差がのっかった状態で説明変数が生成されている」ということを意味していますから、アヒル本でこのモデルをみたときには興奮して、鼻血がでるんじゃないかと思いました。「そんなうまい話、、」というあれもあったので、鼻血はでなかったわけなのですが。


 

仮想データを発生させて遊んでみる


ということでデータの発生から。先に記述したモデルどおりです。

これをlmで分析すると下の結果。真の値 (t) を使って回帰分析をすると、(当たり前ですが) 切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。説明変数を使って分析をすると、切片は3.04、回帰係数は2.33となりました。

 
 
はたして、説明変数にノイズが含まれていることを仮定したモデルだと、真の値で回帰したときの値に近い推定ができるのでしょうか?発生させた仮想データに対して、説明変数にノイズが含まれていることを表現した以下のモデルで推定します。モデルのコードはそのまんまですね。ここでは説明変数のノイズの大きさに、データ発生の際に設定した値である1を固定値として与えていますが、実際に使うときには、先行研究の結果に基づいて分布を与えることになるのではないかと思います。

 

推定結果はこちら。見事にノイズを含む説明変数をつかってlmで分析したときと同じ値が得られており、回帰係数については設定した値である5とはほど遠い値が得られました。やはりそんなうまい話はないのか。。


 
 
ちょっと普通にX1で回帰分析するモデルで推定してみます。回帰分析における誤差の部分 (e) が、説明変数にノイズがあることを仮定するモデルよりも大きくなっています (ただし、lmで分析したときの残差の標準偏差は3.61なので、不当に大きいとか、そういうことではないです)。「説明変数にノイズが含まれている」という仮定は回帰係数の推定にではなく、誤差の推定に効いてくるようです。

 

よろしい、ならば3項目だ


変数が1つだけというのはさすがに無茶だったのかなと思い、変数を増やしてみました。データの発生は変数の数が増えただけで、先ほどと同様です。また、真の値 (t) を使って回帰分析をすると、やはり切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。3つの変数を平均した値を説明変数としてもちいると回帰係数は3.81となりました。


さて、このデータを以下のモデルで推定します。3項目の平均値をつかうモデルでも推定し、比較を行います。


結果はこちら。やはり、両モデルの違いは、切片や回帰係数の値ではなく、誤差の部分にあらわれています。


なお、このデータを普通にSEMで分析してみると、設定した値5に近い値が得られました。

 
 

まとめ


心理学のデータでは説明変数にノイズが含まれていることが多く、「ノイズを除去してより正確な変数間の推定ができれば最高だ!」と興奮したわけですが、回帰係数の推定に効くという感じではなさそうです (仮想データの発生の仕方が悪いのかなぁ)。このあたり、自分はまだまだモデルについての感覚が甘く、遊んでみないとわからないという部分でした。先行研究の因子分析の結果をつかって、「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思ってます。

 
 

【追記】事前分布の設定をするとうまくいく!?


2016.12.12追記
関西学院大学の清水裕士先生に事前分布の設定をするとうまくいくというご助言をいただきました。
モデルにはパラメタの事前分布についての設定を追加しています。特に、パラメタtについては「平均mu、標準偏差sigmaに従う正規分布から発生している」という仮定をいれており、muとsigmaを新たにパラメタとして推定しています。

 
 

推定結果はこちら。データ発生の際に設定した値はa=3, b=5, e=0.5, mu=0, sigma=1なのですが、(特にtについての) 事前分布を設定しない場合に比べ、今回のモデルではそれぞれ近い値が得られています。ただし、今回は反復回数を10万、間引き数を5に設定していて、4万のmcmcサンプルがあるのですが、有効サンプル数 (n_eff) がとても少ないです。もともと説明変数にノイズが含まれていることを表現したモデルだと反復回数や間引き数を多くしないと収束していなかったのですが、事前分布の設定を追加すると、さらに収束のために反復回数や間引き数が必要な感じです。

春の方法論セミナーに参加してきました


日本社会心理学会の第3回春の方法論セミナーに参加してきました。
第3回のテーマは「統計的因果推論への招待」で、過去回と同じく、説明がわかりやすくてとても楽しかった。
以下は感想。これから勉強。

統計的因果推論では、「AがBの原因である」ということの条件に、AがあるときのBの確率(p(B|A))とAがないときのBの確率(p(B|¬A))で、「p(B|A) > p(B|¬A)」ということがある (ようだ)。この条件が必要条件と考えられているのか、十分条件と考えられているのかはわからない。

でも、「p(B|A) > p(B|¬A)」だけだと、「AがBの原因である」場合と「BがAの原因である」場合の区別がつけられない。なので、「p(B|A) > p(B|¬A)」は「AがBの原因である」ことの必要条件であって、十分条件ではない 。ということで、「AがBの原因である」というためには、「p(B|A) > p(B|¬A)」とは別の基準を満たす必要がある。たとえば、時間的先行性で、AがBよりも時間的に先行している必要がある、とか。あとは、AとBの共通原因の影響を排除できているかどうか、とか、だろうか。とりあえず、僕が持ってる素朴な因果概念はこの3つが満たされていたらOKだと思う。

個人的に目から鱗だったのが、星野先生のトークで出てきた「割り当てられなかった条件のデータは欠測データになる」ということ。

大塚先生のトークにあったように、「(Aさんは) 予防接種したからインフルにかからなかった」は「予防接種してインフルにかからなかった」と「予防接種しなくてインフルにかかった」を観察しないといけない。けど、現実ではどちらかしか観察できない。個別事象の因果についてはその検証をするために可能世界を想定しないといけない。心理学のように、個人レベルではなく個体群レベルで考えるなら可能世界を想定する必要はないのかもしれないけど、個体群レベルで考えてもやはり、割り当てられなかった条件のデータは (観察したいけど) 観察できず欠測データになる、というのは同じ (だろう)。

実験で無作為割り付けが重要なのは、実験操作のみが群間で異なるようにするため。だからこそ、従属変数の差を操作に帰属できる。もちろん、これは当たり前のことなんだけど、知識としてもっていただけというかなんというか、まぁ、「無作為割り付けしてるからいいじゃん」みたいな感覚だった。無作為割り付けは交絡要因を統制する消極的な方法、これも知識としてもってる。でも、「本当に実験操作のみが群間で異なるのですか?」といわれたら「せやかて!」と思うレベルには正当化している。

前置きが長くなったけど、要は、「割り当てられなかった条件のデータは欠測データになる」と聞いた時に、観測データを欠測セルに代入している操作が頭に浮かんできて (それ自体は間違ったイメージだろうけど)「うげっ」となった、という話です。傾向スコア万歳。積極的に共変量を除去していこう。

とはいっても、共変量を統計的に除去する、というのには昔から何か引っかかるものがあって、もうちょっと自分の考えを整理しないといけないなぁと思っているところ。でも、たぶん、実際の使用に引っかかってるだけのような気がする。共変量として何を入れるべきかということについて、従属変数と関連が強いが欠測と関連が低い変数を、といった基準もあるのだろうけど、概念的に共変量として考えられるか、ということも重要だろう。媒介変数として機能しているものを共変量として扱っちゃったらまずいんでしょう?「共変量」って言葉の (個人的) イメージが悪いのかもしれない。あと、共変量をめちゃくちゃ除去するという方法は、因果関係があるかないかの判断には使えると思うけど、そこで推定される値に心理学的な意味はあるのだろうか、という疑問。

あと思ったことは、
1) 測定の信頼性についてちゃんと勉強しよう、
2) 林先生のトークであった「理が知りたい」と「介入がしたい」の違いがわかってない、
3) 講演概要にあった「因果推論と統計推論とは、そもそも世界と認識についての考え方 (パラダイム) が根本的に異なっている」というのがわかってない、
の3本。

昨年のGLMMに引き続き、わかりやすく説明してもらったので、関連本を読める気がしてきた。
ありがとうございました。

広島ベイズ塾で発表しました


魁!!広島ベイズ塾でベータ分布について発表しました。今回の会合では「分布感をつかむ」といういまいちつかみきれないテーマで、いくつかの分布についてお勉強しました。正規分布、ポアソン分布、ベータ分布、ガンマ分布、コーシー分布、負の二項分布です。ガンマ分布、コーシー分布については、コワイ本読書会のなかでもよく事前分布としてつかわれていたんですが、どんな分布かまったく把握していなかったのでとても勉強になりました。

負の二項分布についてはあらためて復習しないとわからないですね。「2回くらい変身します」っていわれたときのフリーザ感ったらなかったですね。これが分布感か。

分布の勉強とは離れて、ベイズについてもいろいろ教えてもらいました。なにやら変分ベイズ(?)やノンパラベイズ(?)がrstanで近々できるようになるそうです(ノンパラベイズはもうちょっとかかるんだったかな)。変分ベイズができればmixtureがMCMCできるらしい。ノンパラベイズができればmixtureのときの妥当なクラスの数も推定可能になるらしい。原理も教えてもらったのですが、まったくわかりませんでした。でも、期待で胸が膨らみます。

一応、ベータ分布についての資料をあげておきます。


ハムストリングスと腰痛

ハムストリングスはふとももの裏の筋肉の総称です。大腿二頭筋、半腱様筋、半膜様筋で構成されています。股関節と膝関節をつなぐ二関節筋で、基本的には膝関節屈曲動作、股関節伸展動作ではたらきます。

この柔軟性が低下すると、わたしの場合(おそらく一般的にも)、腰痛がでやすいです。ハムストリングスは股関節についていますので、ここが硬くなると上半身を前傾させるときに、骨盤からうごいてくれなくなります。そうすると、腰の部分を大きく曲げる必要があるため、負担がかかって痛みにつながるということらしいです。

ということで、腰痛が気になる人は、一度ハムストリングスの柔軟性をチェックしてみてはいかがでしょうか。もちろん、腰痛の原因は他にもたくさんありますので、ここの柔軟性にだけ注意すればOKというわけではありませんが。

もも裏のストレッチといえば前屈が一般的ですが、「手をどこまで伸ばせるか」ということを意識しすぎて、意外ともも裏が伸びていないということがよくあります。しっかり骨盤から倒してやって、もも裏を伸ばすことを意識する必要があります。しかしながら、「ハム伸ばさないとなー」と思うころには、ある程度ハムストリングスの柔軟性は低下しており、前屈の際にしっかり骨盤から倒すということがすでに難しくなっている、ということが少なくありません。場合によっては、無理に前屈をすることで腰の負担を高めてしまうかもしれません。

おすすめは、しゃがんだ状態で足首をつかんで膝を伸ばしていく、という方法です。このとき、胸とふとももを密着させておくことに注意して、膝を伸ばしていきます。普段、柔軟体操をしていない方は、膝を少し伸ばそうとするだけで、強烈なストレッチを体感できることと思います。あくまで「ハムストリングスを伸ばす」ことが目的ですので、膝は伸びきらなくてOKです。胸とふとももが離れない程度に膝を伸ばしてください。文字で書くとわかりづらいですね。ジャックナイフストレッチという呼び方があるようなので、詳しくは検索していただければと思います。

立位での前屈との違いですが、「前屈は膝関節を伸ばした状態で股関節を屈曲させる。先ほどの方法は股関節を屈曲させた状態で膝関節を伸ばす」という違いがあります。膝関節を伸ばす際には、大腿四頭筋(ふとももの前側) の力を使えますので、ハムストリングスを伸ばしやすいです。