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

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

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

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に置換

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

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

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

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

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

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

A(isnan(A))=0

としてやるだけです。

これだけの操作で

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

マーカーセットの作成

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

前回の記事では、「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で開いて確認しながらすると良いでしょう。
あとはひたすら根気よく調整を繰り返すのみです。

モデルのダウンロード

今回はモデルの作成方法について解説していきたいと思います。

まずは、サンプルモデルをダウンロードしましょう。

SimTK: OpenSim: Downloads


f:id:simulate:20150523190027p:plain

Models.zipをダウンロードすると、基本的なモデルがほとんどそろっています。

まずは、歩行分析などに用いられるモデルを開いてみましょう。

f:id:simulate:20150523190555p:plain

Fileから「Open Model」を選択します。
そして、先ほどダウンロードした、Models.zipの中の「Gait2354_Simbody」フォルダを開き、「gait2354_simbody.osim」を開きます。

f:id:simulate:20150523190857p:plain

そうすると、このような腕のないヒトのモデルが開きます。

f:id:simulate:20150523190952p:plain


ここにモーションキャプチャで測定したデータを読み込むためのマーカーセットを作成する必要があります。
そのためにデフォルトのマーカーセットをまずは確認しましょう。

f:id:simulate:20150523191146p:plain

f:id:simulate:20150523191227p:plain


同じフォルダ内に「gait2354_Scale_MarkerSet.xlm」というファイルがありますので、これを開きます。
そうすると、このようにマーカーが表示されると思います。

f:id:simulate:20150523191354p:plain

ここから、このマーカーを自分の実験データに合わせて動かしていくわけですが、それは次回にしたいと思います。

OpenSim

OpenSimとは、無料で筋骨格モデルを使ってバイオメカニクス的分析を行えるソフトのことです。

無料であるのは大変ありがたいことですが、残念なことにマニュアルのほとんどが英語で書かれており、使い方を覚えることが困難です。

そこで、ここではその使い方を解説しようと思います。

相関関係の検討

Matlabには標準で簡単な統計を行える関数が備わっています.

最も代表的なものは相関分析でしょう.

corrcoef関数を使うことでPearsonの積率相関係数が算出できます.

その他にもstatistics toolboxを購入することで,t検定などの分析もできるようになります.

statistics toolboxに興味がある方がググってもらうとして,相関分析の方法を説明します.

 

[r,p]=corrcoef(変数)

と書いて使用します.

rは相関係数,pは有意確率です.

変数には,系列を行に入れていきます.例えば,1行目に疾走速度,2行目に腿上げ角度,という具合に.

あとは実行してやると相関行列と有意確率の行列が出てきます.

 

ちなみに,変数は何系列あっても大丈夫です.

私が試した限りでは,30系列くらいは問題なく実行されました.

おそらく,使用しているマシンパワーに依存していると思われるので,相当な数でも分析してくれると考えられます.

Matlab 起動時のフォルダを設定

こんにちは.

今日は起動時に開かれるフォルダの設定方法を紹介したいと思います.

Matlab起動時のフォルダは,Windowsではインストール時に作成されるMydocument内のMATLABフォルダ,Macではユーザーフォルダに設定されていると思います.

これを別のフォルダに設定したい場合,いくつか方法があります.

1.userpath関数

まず紹介するのは,userpathという関数を使う方法です.

例えば,Machogeというユーザーのフォルダの中の「analysis」フォルダにmatlabのデータを置いて使うことを想定した場合には


userpath('/Users/hoge/analysis/') ;

となります.

これをコマンドウインドウに入力して実行した後にMatlabaを再起動してみてください.

 

2.startup.m

Matlabの初期設定の起動フォルダに「startup.m」を作成する方法です.
startup.mというファイルに


cd('/Users/hoge/analysis/') ;

と書いて保存しておくだけです.

この方法では.startup.mに書いておくだけでMatlabのあらゆる初期設定ができます.

 

他にも方法はあるのかもしれませんが,私が知っているのはこの2つだけです.
ぜひ試してみてください.
今後も役立つMatlabの使い方や設定方法などを紹介していきたいと思います.