線形混合効果モデルのデータの準備
テーブルとデータセット配列
線形混合効果モデルを当てはめるには、データをテーブルまたはデータセット配列に保存しなければなりません。テーブルまたはデータセット配列には、応答変数を含む各変数用の列がなければなりません。より具体的には、テーブルまたはデータセット配列 (tbl
) には、次が格納されていなければなりません。
応答変数
y
連続変数またはグループ化変数である予測子変数
Xj
グループ化変数
g1
、g2
、...、gR
ここで Xj
および gr
のグループ化変数は、categorical 配列、logical 配列、文字配列、string 配列、または文字ベクトルの cell 配列です (r = 1、2、...、R)。
各行が観測値を表すようにデータを整理しなければなりません。また、各行には変数の値とその観測値に対応するグループ化変数の水準が含まれていなければなりません。たとえば、4 つの処理オプションを用いる実験データがある場合、個体の母集団 (ブロック) から無作為に選択した 5 つの異なる種類の個体のテーブルまたはデータセット配列は、次のようになります。
ブロック | 処理 | 応答 |
---|---|---|
1 | 1 | y11 |
1 | 2 | y12 |
1 | 3 | y13 |
1 | 4 | y14 |
... | ... | ... |
5 | 1 | y51 |
5 | 2 | y52 |
5 | 3 | y53 |
5 | 4 | y54 |
トマトの苗木を異なる 4 種類の肥料で栽培する効果を調査する分割実験についてここで検討します。トマトの苗木が栽培されている土壌はその種類に基づいて 3 つのブロックに分割されています。砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。次に、プロット内でトマトの苗木はサブプロットに分割され、それぞれのサブプロットに 4 つの肥料のいずれかが処理されます。この実験によるデータは次のようになります。
土壌 | トマト | 肥料 | 収穫数 |
---|---|---|---|
'砂質' | 'プラム' | 1 | 104 |
'砂質' | 'プラム' | 2 | 136 |
'砂質' | 'プラム' | 3 | 158 |
'砂質' | 'プラム' | 4 | 174 |
'砂質' | 'チェリー' | 1 | 57 |
'砂質' | 'チェリー' | 2 | 86 |
... | ... | ... | ... |
'砂質' | '枝付き' | 3 | 99 |
'砂質' | '枝付き' | 4 | 117 |
'シルト' | 'プラム' | 1 | 120 |
'シルト' | 'プラム' | 2 | 115 |
... | ... | ... | ... |
'粘土質' | '枝付き' | 3 | 111 |
'粘土質' | '枝付き' | 4 | 105 |
fitlme
の入力引数 formula
を使用して近似させるモデルを指定しなければなりません。
一般に、モデル仕様の式は 'y ~ terms'
という形式の文字ベクトルまたは string スカラーです。線形混合効果モデルでは、この式は 'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)'
の形式になります。ここで fixed
には固定効果の項が、random1, ..., randomR
には変量効果の項が含まれます。たとえば前の肥料の実験では、次の混合効果モデルについて検討します。
ここで i = 1、2、...、60 で、インデックス m は肥料の種類、j はトマトの種類に対応し、k = 1、2、3 はブロック (土壌) に対応します。Sk は k 番目の土壌の種類を表し、I[F]im は肥料の水準 m を表すダミー変数です。同様に、I[T]ij はトマトの種類の水準 j を表すダミー変数です。
式 'Yield ~ 1 + Fertilizer + Tomato + (1|Soil)+(1|Soil:Tomato)'
を使用してこのモデルを当てはめることができます。
式を使用したモデルの指定方法の詳細については、式と計画行列の関係を参照してください。
計画行列
式を使用したモデルの記述が困難な場合は、固定効果と変量効果を定義する計画行列を作成し、fitlmematrix(X,y,Z,G)
を使用してモデルを当てはめることができます。計画行列は次のように作成しなければなりません。
固定効果と変量効果の計画行列 X
および Z
の作成方法は次のとおりです。
ones(n,1)
を使用して切片に 1 の列を入力します。ここで n は観測の総数です。X1
が連続変数の場合は、別の列にそのままX1
を入力します。X1
が m 水準のカテゴリカル変数の場合、X
内においてX1
の m – 1 水準には m – 1 個のダミー変数が存在しなければなりません。たとえば、4 つの異なる供給業者からの原材料の品質が、生産ラインの生産性に与える影響について調査する実験を考えます。切片と供給業者を固定効果の項とする線形混合効果モデルを当てはめ、切片は変量効果の項とし、基準コントラスト コーディングを使用する場合、固定効果および変量効果の計画行列を次のように作成しなければなりません。
D = dummyvar(provider); % Create dummy variables X = [ones(n,1) D(:,2) D(:,3) D(:,4)]; Z = [ones(n,1)];
基準コントラスト コーディングは最初の供給業者を基準として使用します。モデルには切片があるため、残りの 3 つの供給業者だけにダミー変数を使用しなければなりません。
予測子変数
X1
とX2
の交互作用項が存在する場合は、ベクトルX1
とX2
の要素ごとの積で形成される列を入力しなければなりません。たとえば、長期間の調査において、切片、連続処理因子、連続時間因子およびそれらの交互作用が固定効果として含まれており、時間が変量効果の項であるモデルを当てはめる場合は、固定効果および変量効果の計画行列は次のようになります。
X = [ones(n,1),treatment,time,treatment.*time]; y = response; Z = [time];
グループ化変数 G
の作成方法は次のとおりです。
グループ化変数ごとの 1 つの列と、入れ子の場合はグループ化変数の要素ごとの積の列があります。
たとえば、ブロック (block
) 内のプロット (plot
) をグループ化する場合は、block
別に plot
の要素ごとの積の列を追加しなければなりません。より具体的には、ブロック分割試験において切片および連続処理因子が固定効果として含まれるモデルを当てはめる場合、切片と処理がブロック内で入れ子にされたプロットによってグループ化されていると、計画行列は次のようになります。
X = [ones(n,1),treatment]; y = response; Z = [ones(n,1),treatment]; G = [block.*plot];
前述した原材料の品質の例では、原材料がバルクで届き、バルクは供給業者の入れ子になると仮定します。線形混合効果モデルを当てはめる場合、切片が供給業者内のバルクによってグループ化されるのであれば、計画行列は次のようになります。
D = dummyvar(provider); X = [ones(n,1) D(:,2) D(:,3) D(:,4)]; y = response; Z = ones(n,1); G = [provider.*bulks];
前の長期間の調査の例で、調査に参加する被験者によってグループ化された切片と時間の変量効果を追加する場合は、計画行列は次のようになります。
X = [ones(n,1),treatment,time, treatment.*time]; y = response; Z = [ones(n,1),time]; G = subject;
行列形式とテーブルおよびデータセット配列の関係
次に示すように fitlme(tbl,formula)
と fitlmematrix(X,y,Z,G)
の機能は同じです。
y
は n 行 1 列の応答ベクトルです。X
は n 行 p 列の固定効果計画行列です。fitlme
はこれをformula
内の式fixed
から作成します。Z
は R 行 1 列の cell 配列です。Z{r}
は n 行 q(r) 列の変量効果計画行列であり、formula
内のrandom
の r 番目の式から作成されます (r = 1、2、...、R)。G
は R 行 1 列の cell 配列です。G{r}
は n 行 1 列のグループ化変数g
r であり、formula
内の M(r) 水準またはグループにあります。
たとえば、tbl
が応答変数 y
、連続変数 X1
および X2
、グループ化変数 g
を含むテーブルまたはデータセット配列の場合に、fitlmematrix(X,y,Z,G)
を使用して式 'y ~ X1+ X2+ (X1*X2|g)'
に対応する線形混合効果モデルを当てはめるには、その入力引数を次のように対応させなければなりません。
y = tbl.y X = [ones(n,1), tbl.X1, tbl.X2] Z = [ones(n,1), tbl.X1, tbl.X2, tbl.X1.*tbl.X2] G = tbl.g
参考
LinearMixedModel
| fitlme
| fitlmematrix