線形混合効果モデルのデータの準備

テーブルとデータセット配列

線形混合効果モデルに近似するには、データをテーブルまたはデータセット配列に保存しなければなりません。テーブルまたはデータセット配列には、応答変数を含む各変数用の列がなければなりません。より具体的には、テーブルまたはデータセット配列 (tbl) には、次が格納されていなければなりません。

応答変数y
連続変数またはグループ化変数である予測子変数X_j
グループ化変数g₁、g₂、...、g_R

ここでX_jおよびg_rのグループ化変数は、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には変量効果の項が含まれます。たとえば前の肥料の実験では、次の混合効果モデルについて検討します。

$y_{i m j k} = β_{0} + \sum_{m = 2}^{4} β_{1 m} I {[F]}_{i m} + \sum_{j = 2}^{5} β_{2 j} I {[T]}_{i j} + b_{0 k} S_{k} + b_{0 j k} {(S * T)}_{j k} + ε_{i m j k},$

ここで i = 1、2、...、60 で、インデックス m は肥料の種類、j はトマトの種類に対応し、k = 1、2、3 はブロック (土壌) に対応します。S_kは 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 variablesX = [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