読者です 読者をやめる 読者になる 読者になる

スポーツバイオメカニクス MatlabとOpenSIM

スポーツバイオメカニクスの分析によく用いられるMatlabとOpenSIMの使用に関する備忘録

バイオメカニクスの教科書

今日は、スポーツバイオメカニクスを勉強するにあたって、おすすめの教科書を紹介します。

まずはこちら。

阿江通良・藤井範久 著
「スポーツバイオメカニクス20講」


大学でバイオメカニクスを学んだ人で知らない人はおそらくいないと思います。
幅広い内容が簡潔に、たっとこれだけのサイズにまとめられている素晴らしい教科書です。
ただし、初心者がこの本で本格的に勉強するには、詳細な説明が少なすぎてやや難しいかもしれません。


次に紹介するのは、
ウインター 著
バイオメカニクス 人体運動の力学と制御」

こちらも非常に有名な教科書です。
もともとは英語で書かれた本ですが、日本語に訳されました。
初版はかなり古い本ですが、こちらは原著第4版ということで、モーションキャプチャなど、最近の新しい情報も大幅に書き加えられています。
また、信号処理についても解説があり、これからバイオメカニクスを勉強する方にとって非常に良い本だと思います。


3冊目は
ゴードン 著
「身体運動のバイオメカニクス研究法」

こちらは、最初に紹介した阿江通良先生が訳した本です。
この本もバイオメカニクスの教科書として、日本ではスタンダードな1冊です。
ウインターの本よりも価格が安いので、学生にも購入し易いのは魅力的です。

他にもおすすめの教科書はいくつかありますが、今日はこの3冊の紹介に留めておきたいと思います。

Matlab Academy —Matlabの使い方を基礎から—

Matlab

Matlabの基礎的な使い方はどうやって勉強しましたか?

本を読んだり、ネットで調べたりして勉強したという人が多いと思います。
ネットで調べると、かなりわかりやすいサイトも増えてきましたが、基礎の基礎を学ぶのはなかなか大変です。
おそらく、今、何をやっているのかが理解しにくい作業だからでしょう。

これからMatlabを始めようという方に、強い味方ができました。

MathWorksが提供する「Matlab Academy」です。

これまでは英語でしか提供されていなかったのですが、日本語でも利用できるようになりました。

高度なプログラムは有料提供となっていますが、入門コースは無料で提供されています。
ブラウザ上で、Matlabが動くので、Matlabを持っていなくても利用できます。
使用するためには、MathWorksアカウントに登録するだけです。

入門コースをマスターするだけで、Matlabに関する情報の多くが理解でき、利用できるようになります。
これからMatlabを使えるようになりたいという方はぜひ利用してみてください。

t検定

Matlab

毎度のことながら、久々の更新です。

以前、分散分析については記事にしていますが、分散分析よりももっとシンプルで用いられることも多いt検定について書いていませんでした。
そこで、今回はMatlabでt検定を行う方法を紹介したいと思います。

対応ありt検定(1標本)

そのまんまですが、対応ありの場合に用いるt検定です。
同じ被験者が異なる2つの条件で試技を行った際に比較する方法です。
条件1のときの値をa、条件2のときの値をbとしたとき、

h=ttest(a,b);

となります。
有意水準は5%に設定されており、有意差が認められた場合には、h=1と出力され、有意差が認められなかった場合には、h=0と出力されます。

有意確率も確認したい場合には、

[h,p]=ttest(a,b);

とすることで、pに有意確率が出力されます。

また、論文を書く際には、t値も必要になるかと思いますが、t値も含めた詳細な統計量を出力する際には

[h,p,ci,stats]=ttest(a,b);

とすることで確認できます。

MathWorksの公式な説明はこちらにあります。

対応なしt検定(2標本)

異なる2群間の差の検定には、対応なしのt検定を用います。
こちらも方法は至ってシンプルで、

[h,p,ci,stats]=ttest2(a,b);

とするのみです。
対応あり要因は「ttest」関数を用いましたが、対応なし要因は「ttest2」を用います。
利用方法は基本的に同じですので、割愛します。
また、対応なしt検定のMathWorksの説明はこちらから参照できます。

最後に

差の検定を行うにあたり、「とりあえず」で簡単に行えるのは今回紹介したt検定です。
MatlabSPSSを使わずとも、Excelの関数でも簡単に行えます。
ただし、検定にかけたい変数が膨大な場合には、Matlabで検定を行うことが便利な場合が多いでしょう。
また、Matlabを用いる利点としては、SPSSと比較して有効数字が大きいことが挙げられます。
特に、小数点以下にゼロが何桁も続くような非常に小さい値を用いる場合には、SPSSでは四捨五入されてしまうことがありますが、Matlabを用いた場合は浮動小数点計算が行われますので、かなりの精度で計算してくれます。
これは、Matlabを統計に用いることの大きな利点と言えるでしょう。

ベンチマーク

Matlab

今日はMatlabベンチマークについて紹介したいと思います。

ベンチマークと言えば「Geekbench」などが有名ですね。

Geekbench

Geekbenchは、いろいろなプラットフォーム(OS)で使えるのが最大の利点です。
私もPCを購入する際のひとつの観点として、Geekbenchのスコアは気にしています。

しかし、Matlabで計算をする人にとって、そんなベンチマークのスコアよりもMatlabでの計算時間がどれくらい変わるのかの方が気になるのではないでしょうか。

そんなときにおすすめしたいのがMatlabベンチマーク機能です。
使い方は、コマンドウインドウに

bench
と入力するだけです。

 
このベンチマークを使用する際の注意点としては、異なるバージョンでは計算方法が異なるらしく、直接結果を比較することができないそうです。

バージョンが異なると比較できないとなると、参考になるのか微妙な点もありますが、参考にどうぞ。

反復測定2元配置分散分析 【Statistics and Machine Learning Toolbox】

Matlab

Statistics and Machine Learning Toolbox

MatlabのStatistics Toolboxがバージョンアップされて、「Statistics and Machine Learning Toolbox」になりました。
これまでは反復測定分散分析はできなかったようですが、このバージョンアップで正式に対応したようです。

これまでのStatistics Toolboxでは、二元配置分散分析を行うことはできたのですが、反復のない分散分析のみ対応していました。
ここで、「反復」って何だ?と思う方もいるでしょう。
詳細は統計の本やサイトに譲りますが、「反復測定」というのは、いわゆる「対応あり」のことを言います。
つまり、これまでのToolboxでは対応なしの分散分析しかできなかったわけです(多分)。

この度新しく追加された反復測定分散分析ですが、これまでの反復のない分散分析と比較して、使い方が非常に難しいのです。

反復測定分散分析 —MATLAB—
正直なところ、Mathworksの解説を読んでもさっぱりわかりません。
私は統計の専門家ではないので、詳しいことはわかりませんが、説明文の登場する用語たちも決して統計の専門用語ではないようで、わからない言葉を検索してもMathworksのページばかりが登場します。

そこで、今回はMatlabを用いた反復測定分散分析の使い方の解説をしようと思います。

データセットの準備

今回はランダムな数字を発生させてデータを用意しました。

ある試技をAとBの2条件で同じ被験者に行わせたとします。
そして、時刻1と時刻2のときのあるパラメータを算出しました。

したがって、同じ被験者から得られたデータを用いて、条件(A, B)×時刻(1, 2)の2要因の分散分析を行います。

データはSPSSで分析する際と同じように、1行ごとに1人の被験者のデータを入力します。
列は条件になります。

f:id:simulate:20170131195111p:plain

これをもとに、データセットのテーブルを作成します。

T=table(data(:,1),data(:,2),data(:,3),data(:,4),'VariableNames',{'A_1','A_2','B_1','B_2'});

そうすると、このようにテーブルができます。

f:id:simulate:20170131195157p:plain

条件の指定

データセットが出来たら、今度は条件を指定したテーブルを作成してやります。
今回は条件A、Bと時刻1、2です。

まずは、セル配列として、各列がどの条件にあてはまるのかを指定します。
ここでは、データセットの列は、セル配列の行にあたります。

w1={'A','A','B','B'}';w2={'t1','t2','t1','t2'}';

 2要因ですので、2つの指定が必要です。
1つめの要因では、条件AとBを指定しています。
2つめの要因では、時刻1と時刻2を指定しています。

つまり、データセットテーブルの1列目は条件Aの時刻1、2列目は条件Aの時刻2、3列名は条件Bの時刻1、4列目は条件Bの時刻2のデータであることを示しています。

これをさきほどの容量でテーブルにしてやります。

within=table(w1,w2,'VariableNames',{'Condition','Time'});

反復測定モデルの近似

ここまでできたら、分散分析の内容に入っていきます。
anovaにかける前に、反復測定モデルを作成しなければなりません。
Mathworksのページでは、このような解説がされています。

反復測定モデルの近似 —MATLAB—

では、具体的に説明していきます。
必要なのはたった1行のプログラムです。

rm=fitrm(T,'A_1-B_2~1','WithinDesign',within);

引数を説明すると、
「T」はデータセットテーブルです。
「'A_1-B_2~1'」はモデル仕様の式です。
これに関して、Mathworksの解説を読んでもさっぱりわかりません。
ウィルキンソンの〜と解説されていますが、何のことやら。

モデル仕様は「~」の左側と右側に分けて考えます。
左側が応答、右側が項です。
それでもよくわかりませんね。

基本的に、反復測定2元配置分散分析では、項は「1」でOKです。
応答には、データセットの中の必要な変数をテーブルの変数名で指定してやります。
「-」は範囲を指定する書き方です。
今回の場合、「A_1-B_2」と表記しています。
データセットには、左から順にA_1、A_2、B_1、B_2と入っていますから、この書き方をすると「A_1」から「B_2」までの変数を使うことになるわけです。

A_1とB_2だけを使う、というような場合には、「A_1,B_2」と書けばいいようです。
とはいえ、面倒なので分散分析のデータセットには不要なデータは入れずに、「左端の変数名-右端の変数名」で表記できるようにしましょう。

次に、「'WithinDesign'」と「within」ですが、被験者内要因の定義を示します。
「'WithinDesign'」で反復測定モデル(つまり、被験者内要因があること)であることを宣言し、その内容を次の引数として入力します。
「within」は条件の指定のところで作成したテーブル名です。

ここまで行うと、「rm」という変数ができます。

f:id:simulate:20170131195248p:plain

これは、この先使っていきますが、中身を見る必要はいまのところありません。
「ああ、できてるな」と思っていただければ結構です。

ANOVAの実行

ついに、分散分析を実行します。
対応なしの二元配置分散分析では、「anova2」を用いますが、対応ありの二元配置分散分析では、「ranova」を用います。

[ranovatbl]=ranova(rm,'WithinModel','Condition*Time');

まずは、入力引数から。

「rm」はさきほど作成した反復測定モデルです。
そのまま入れてやれば良いだけです。

「'WithinModel'」は、このあとのどの条件の主効果、交互作用を見るかを指定するためのものです。
「'Condition*Time'」がCondition要因とTime要因について検討することを示しています。
この指定を忘れると表示される結果が全くことなるものになってしまいます。
この「Condition」と「Time」は条件の指定で作成したテーブルの行タイトル(変数名)で指定してやります。

ちなみに、「[ranova]=ranova(rm)」だけで実行した場合の結果はこのようになります。
f:id:simulate:20170131195317p:plain

これでは交互作用も主効果もわかりません。

こちらが正しく指定したときの出力引数「ranovatbl」はです。
分散分析を実行した結果が入っています。

f:id:simulate:20170131195340p:plain

Fの列にF値が、pValueの列にp値が入っています。
7行目には「交互作用」の結果が、3行目には「Condition」の主効果が、5行目には「Time」の主効果が示されています。

pValueの列の右にある3列は球面性検定の結果に基づいて、補正をした場合のp値が示されています。
このあたりの詳細は統計書に譲ります。
なお、球面性検定もMatlabでできます。

今回はその結果を無視して、球面性の仮定ができた場合について見ていきます。

多重比較(交互作用が有意だった場合)

今回の結果では、交互作用が有意になっています。
交互作用が有意であった場合、多重比較の際にもうひとつの要因で補正してやる必要があります。

そのため、多重比較は次のように行います。

mult1=multcompare(rm,'Condition','By','Time');
mult2=multcompare(rm,'Time','By','Condition');

すると、結果は次のようになります。

f:id:simulate:20170131195410p:plain
f:id:simulate:20170131195430p:plain

mult1について見てみると、
1行目は時刻1のときの条件AとBとの間の差を検討しており、p値が0.1758となっていることから、ここには有意差が認められなかったことを示しています。
その他の結果についても同様に見ていきます。

多重比較(交互作用が有意ではない場合)

交互作用が有意ではない場合は、次のように表記します。

mult1=multcompare(rm,'Condition');
mult2=multcompare(rm,'Time');

そして、結果はこのようになります。

f:id:simulate:20170131195458p:plain

f:id:simulate:20170131195533p:plain

交互作用が有意でない場合、条件ごとまとめて比較して良いので、「mult1」では、時刻を問わず、条件AとBとの間に差があるのかどうかを検討しています。
今回の結果では、p値は0.9770であり、有意差が認められなかったことを示しています。

最後に

今回は、Matlabを用いた反復測定分散分析について紹介してきました。
これらの分析はSPSSでも行えますが、シンタックスを用いて単純主効果検定と多重比較について書き足す必要があります。
また、SPSSでは要因が99個を超えると実行できません。

つまり、SPSSを使うと割と面倒なことになります。
これらの分析を1回だけ行うのであれば、大した問題ではありませんが、何十回も行うのであれば、Matlab上でうまく回す方が賢明でしょう。

ただし、最初の1回はSPSS等の統計ソフトの結果と一致するか確認しておくことをお勧めします。

Nanを0に置換

Matlab

いろいろと複雑な計算をやっていると、「Nan」が値として出力されてしまうことがあります。

その原因も明らかでそれを避けられない場合、これを0に置換した方が後の計算に都合が良い場合も多々あります。

そんなときの置換する方法を紹介します。

A=[1 2 Nan; Nan 8 Nan; 9 3 1]

という状態だったと仮定します。

このAの中のNanを一発で0に置換するには、

A(isnan(A))=0

としてやるだけです。

これだけの操作で

Nanはすべて0になっているはずです。

マーカーセットの作成

OpenSIM

前回の更新から随分時間が経ってしまいましたが、今回はマーカーセットの作成について紹介したいと思います。
作成と言っても、もともとあるモデルを編集して作成していきます。

前回の記事では、「gait2354_Scale_MarkerSet.xlm」を開いたところで終わっていたと思います。
今回はこのマーカーセットをもとに編集していきます。

まず、マーカーセットを編集できる状態にしなければならないのですが、マーカーセットの編集にはOpenSimは使いません。
使用するソフトはExcelテキストエディタです。

Excelを起動し、ファイルを開くメニューから「gait2354_Scale_MarkerSet.xlm」を開きます。

f:id:simulate:20150823175218p:plain


そうすると、このような画面が表示されるはずです。
f:id:simulate:20150823180420p:plain


同じファイルをテキストエディタでも開きます。
(ちなみに、私が使用しているのはサクラエディタです。)
f:id:simulate:20150823180652p:plain

なぜExcelテキストエディタを使用するかというと、Excelは内容を理解しやすいのですが、xmlファイルで正しく保存することができません。そこで、Excelで開いて、内容を理解した上でテキストエディタで編集していこうというわけです。

両方の画面を照らし合わせて見て頂くとわかると思いますが、簡単に解説します。
例えば、「Sternum」のマーカーについて見ていきます。
Excelの画面で「name2」というところにSternumと書かれていると思います。
これは、テキストエディタでは、5行目のという部分にあたります。

同様に、他の項目についても見ていくと、このマーカーはtorsoセグメントを構成していますので、Excelでは「body」の列に「torso」と書かれています。一方、テキストエディタでは、「 torso 」というようにダグ内に書かれています。
それ以降も同様に、Excel画面の1行目の文字が、テキストエディタでは「<>」で囲まれたような状態で書かれていると考えてください。

このようにして、テキストエディタの内容が理解できるようになったところで、テキストエディタで開いたファイルを名前をつけて保存しておきましょう。
万が一、データが壊れてしまった場合に、もともとのマーカーセットのデータまで壊れてしまわないためです。
(再度、ダウンロードするという方法もありますが、面倒なので)

f:id:simulate:20150823181514p:plain



そこまでできたら、さっそく編集に入っていきます。
編集するにあたり、一番重要なのは、マーカーの名前です。
つまり、「Marker name="マーカーの名前"」です。
ここに書かれているマーカーの名称とViconデータ内のマーカーの名称が一致している必要があります。

つまり、Vicon Nexusの赤枠で囲んだ部分に表示されている名前と、ここに入力する名前が完全に一致していなければなりません。

f:id:simulate:20150823182001p:plain

この画面のデータはPlugingaitを使っているので、ネット上にマーカーセットのデータが落ちていますが、独自のマーカーセットを使っている場合には、近いマーカーセットのサンプルをもとに自分でそれを編集してマーカーセットデータを作成しなければなりません。

例えば、マーカーを増やしたい場合には、さきほどのテキストエディタの画面でいうところの5行目から33行目までをコピーして、内容を書きかえてやれば増やすことができ、逆に減らす場合には同様に「Marker」で囲まれた部分を削除してあげればOKです。


次に、実際の編集の仕方ですが、基本的には、どのセグメントを構成するマーカーなのかと、そのセグメントにおける位置を決めてやれば良いだけです。
もしかすると、もっと細かいきまりがあるかもしれませんが、まだそこまで理解できていません。既存のマーカーセットにマーカーを少し加える、もしくは少し位置を変えるくらいの編集であればこの程度で問題ないはずです。


まずは、構成するセグメントを決めます。
bodyダグの中にセグメント名を書き込んでください。

マーカーの位置はセグメント座標系の原点からの位置になりますので、セグメントを間違えるとおかしなことになります。
位置の編集は、OpenSimで開いて確認しながらすると良いでしょう。
あとはひたすら根気よく調整を繰り返すのみです。