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

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

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

  • 応答変数 y

  • 連続変数またはグループ化変数である予測子変数 Xj

  • グループ化変数 g1、g2、...、gR

ここで Xj および gr のグループ化変数は、categorical 配列、logical 配列、文字配列、string 配列、または文字ベクトルの cell 配列です (r = 1、2、...、R)。

各行が観測値を表すようにデータを整理しなければなりません。また、各行には変数の値とその観測値に対応するグループ化変数の水準が含まれていなければなりません。たとえば、4 つの処理オプションを用いる実験データがある場合、個体の母集団 (ブロック) から無作為に選択した 5 つの異なる種類の個体のテーブルは、次のようになります。

ブロック処理応答
11y11
12y12
13y13
14y14
.........
51y51
52y52
53y53
54y54

トマトの苗木を異なる 4 種類の肥料で栽培する効果を調査する分割実験についてここで検討します。トマトの苗木が栽培されている土壌はその種類に基づいて 3 つのブロックに分割されています。砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。次に、プロット内でトマトの苗木はサブプロットに分割され、それぞれのサブプロットに 4 つの肥料のいずれかが処理されます。この実験によるデータは次のようになります。

土壌トマト肥料収穫数
'砂質''プラム'1104
'砂質''プラム'2136
'砂質''プラム'3158
'砂質''プラム'4174
'砂質''チェリー'157
'砂質''チェリー'286
............
'砂質''枝付き'399
'砂質''枝付き'4117
'シルト''プラム'1120
'シルト''プラム'2115
............
'粘土質''枝付き'3111
'粘土質''枝付き'4105

fitlme の入力引数 formula を使用して近似させるモデルを指定しなければなりません。

一般に、モデル仕様の式は 'y ~ terms' という形式の文字ベクトルまたは string スカラーです。線形混合効果モデルでは、この式は 'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)' の形式になります。ここで fixed には固定効果の項が、random1, ..., randomR には変量効果の項が含まれます。たとえば前の肥料の実験では、次の混合効果モデルについて検討します。

yimjk=β0+∑m=24β1mI[F]im+∑j=25β2jI[T]ij+b0kSk+b0jk(S*T)jk+εimjk,

ここで 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 列のグループ化変数 gr であり、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

参考

| |

トピック