nlmefit

非線形混合効果の推定

構文

beta = nlmefit(X,y,group,V,fun,beta0)
[beta,PSI] = nlmefit(X,y,group,V,fun,beta0)
[beta,PSI,stats] = nlmefit(X,y,group,V,fun,beta0)
[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0)
[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0,'Name',value)


説明

beta = nlmefit(X,y,group,V,fun,beta0) は、非線形混合効果の回帰モデルを当てはめ、beta に固定効果の推定を返します。既定の設定では、nlmefit は、各パラメーターが固定効果と変量効果の和となり、変量効果が無相関となる (変量効果の共分散行列が対角行列である) モデルを当てはめます。

X は、h 予測子における n 個の観測の n 行 h 列の行列です。

y は、n 行 1 列の応答ベクトルです。

group は、観測値の m 個のグループを示すグループ化変数です。group は、カテゴリカル変数、数値ベクトル、各行がグループ名に対応する文字行列、string 配列、または文字ベクトルの cell 配列です。グループ化変数の詳細は、グループ化変数を参照してください。

V は、g 個のグループ固有の予測子の m 行 g 列の行列または cell 配列です。これらの予測子は、グループ内のすべての観測に対して同じ値を取ります。V の行は、grp2idx(group) で指定された順序に従って、grp2idx を使用して、グループに割り当てられます。グループ予測子のサイズがグループ全体で変化する場合、V に cell 配列を使用します。グループ固有の予測子がない場合は、V に [] を使用します。

fun は、予測子値とモデル パラメーターを受け入れ、近似値を返す関数のハンドルです。fun は、以下の形式になります。

yfit = modelfun(PHI,XFUN,VFUN)

引数は、次のようになります。

• PHI — モデル パラメーターの 1 行 p 列のベクトル。

• XFUN — 予測子の k 行 h 列の配列。ここで、k は次のようになります。

  • k = 1 (XFUN が X の単一行の場合)

  • k = n_i (XFUN がサイズ n_i の単一グループに対する X の行を含む場合)

  • k = n (XFUN が X のすべての行を含む場合)

• VFUN — グループ固有の予測子。次のいずれかで与えられます。

  • 単一グループまたは V の単一行に対応する 1 行 g 列のベクトル。

  • n 行 g 列の配列。ここで、次の場合、j 行目は V(I,:) です。j 番目の観測がグループ I に含まれる場合。

  V が空の場合、nlmefit は 2 つの入力で modelfun を呼び出します。

• yfit — 近似値の k 行 1 列のベクトル

PHI または VFUN のいずれかが単一行を含む場合、その行は他の 2 つの入力引数のすべての行に対応します。

メモ

modelfun が 1 回の呼び出しでモデル パラメーターの複数のベクトルに対して yfit を計算できる場合、'Vectorization' パラメーター (後述) を使用するとパフォーマンスが向上します。

beta0 は、q の固定効果に対する初期推定をもつ q 行 1 列のベクトルです。既定の設定では、q はモデル パラメーター数 p です。

nlmefit は、統合した変量効果で周辺尤度への近似を最大化することでモデルを当てはめます。ここで、以下の内容を仮定します。

• 変量効果は多変量正規分布で、グループ間で独立している。

• 観測誤差は互いに独立で同一の正規分布に従い、変量効果から独立している。

[beta,PSI] = nlmefit(X,y,group,V,fun,beta0) は、変量効果に対する r 行 r 列の推定された共分散行列、PSI も返します。既定の設定では、r はモデル パラメーター数 p と等しくなります。

[beta,PSI,stats] = nlmefit(X,y,group,V,fun,beta0) は以下のフィールドをもつ stats 構造体も返します。

• dfe — モデルに対する誤差自由度

• logl — 近似したモデルに対する最大対数尤度

• rmse — 推定誤差分散の平方根 (exponential 誤差モデルの対数スケールで計算)

• errorparam — 誤差分散モデルの推定パラメーター

• aic — 赤池情報量基準。aic = -2 * logl + 2 * numParam により計算されます。ここで、numParam は近似パラメーター (変量効果の共分散行列の自由度、固定効果の個数、誤差モデルのパラメーターの個数など) の個数、logl は構造体 stats のフィールドです。

• bic — ベイズ情報量基準。bic = –2*logl + log(M) * numParam のように計算します。

  • M は、グループの数です。

  • numParam と logl は、aic で定義されます。

  一部の文献では、bic の計算を bic = -2*logl + log(N) * numParam としています。N は観測数です。

• covb — パラメーター推定の推定された共分散行列を返します。

• sebeta — 以下に対する標準誤差 beta

• ires — 母集団の残差 (y-y_population)。ここで y_population は、個別の予測値です。

• pres — 母集団の残差 (y-y_population)。ここで y_population は、母集団の予測値です。

• iwres — 個別の重み付きの残差

• pwres — 母集団の重み付きの残差

• cwres — 条件付き重み付きの残差

[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0) は、m 個のグループに対して推定された変量効果の r 行 m 列の行列、B も返します。既定の設定では、r はモデル パラメーター数 p と等しくなります。

[beta,PSI,stats,B] = nlmefit(X,y,group,V,fun,beta0,'Name',value) は、オプション パラメーターの名前と値のペアを 1 つ以上指定します。一重引用符で囲んで Name を指定します。

既定とは異なるモデルを当てはめるには、以下のパラメーターを使用します (既定のモデルは、FEConstDesign と REConstDesign の両方を eye(p) に設定するか、FEParamsSelect と REParamsSelect の両方を 1:p に設定することにより得られます)。'FE' の接頭辞をもつパラメーターと、'RE' の接頭辞をもつパラメーターを最大で 1 つ使用します。関数 nlmefit では、少なくとも 1 つの固定効果と 1 つの変量効果を指定する必要があります。

パラメーター	値
FEParamsSelect	1 ～ p のインデックスの数値ベクトルまたは 1 行 p 列の logical ベクトルとして与えられた、固定効果を含むパラメーターのベクトル PHI の要素を指定するベクトル。q が指定した要素数である場合、モデルは q 個の固定効果を含みます。
FEConstDesign	p 行 q 列の計画行列 ADESIGN。ここで、ADESIGN*beta は PHI の p 要素の固定要素です。
FEGroupDesign	m 個のグループごとに異なる p 行 q 列の固定効果の計画行列を指定する p x q x m 配列。
FEObsDesign	n 個の観測ごとに異なる p 行 q 列の固定効果の計画行列を指定する p x q x n 配列。
REParamsSelect	1 ～ p のインデックスの数値ベクトルまたは 1 行 p 列の logical ベクトルとして与えられた、変量効果を含むパラメーターのベクトル PHI の要素を指定するベクトル。モデルは r 個の変量効果を含みます。ここで、r は、指定した要素数です。
REConstDesign	p 行 r 列の計画行列 BDESIGN。ここで、BDESIGN*B は PHI の p 要素のランダム要素です。
REGroupDesign	m 個のグループごとに異なる p 行 r 列の変量効果の計画行列を指定する p x r x m 配列。
REObsDesign	n 個の観測ごとに異なる p 行 r 列の変量効果の計画行列を指定する p x r x n 配列。

尤度を最大にするための反復アルゴリズムを制御するには、以下のパラメーターを使用します。

Parameter	値
RefineBeta0	nlmefit が変量効果のない modelfun を近似した後、beta0 を beta で置き換えることで、beta0 の初期調整を行うかどうかを指定します。選択肢は、'on' と 'off' です。既定値は 'on' です。
ErrorModel	誤差項の形式を指定する文字ベクトルまたは string スカラー。既定値は 'constant' です。各モデルは、標準正規 (ガウス) 変数 e、関数値 f、1 つまたは 2 つのパラメーター a と b を使用して、誤差を定義します。選択肢は、以下のとおりです。 'constant': y = f + a*e 'proportional': y = f + b*f*e 'combined': y = f + (a+b*f)*e 'exponential': y = f*exp(a*e)、または等価な log(y) = log(f) + a*e このパラメーターが指定された場合、出力される stats.errorparam フィールドの値は以下のようになります。 a ('constant' および 'exponential' の場合) b (次の場合) 'proportional' [a b] (次の場合) 'combined'
ApproximationType	モデルの尤度を近似するために使用する手法。選択肢は、以下のとおりです。 'LME' — beta と B の現在の条件推定で線形混合効果モデルに対する尤度を使用します。これは既定の設定です。 'RELME' — beta と B の現在の条件推定で線形混合効果モデルに対する制限付き尤度を使用します。 'FO' — 変量効果のない 1 次ラプラシアン近似。 'FOCE' — B の条件推定での 1 次ラプラシアン近似。
Vectorization	modelfun への PHI、XFUN、VFUN の入力引数の許容可能なサイズを指定します。選択肢は、以下のとおりです。 'SinglePhi' — modelfun は一度に 1 組のモデル パラメーターしか受け入れることができないので、各呼び出しでは PHI が単一の行ベクトルになっていなければなりません。必要な場合、nlmefit はループ内で modelfun を呼び出しますが、1 回の呼び出しでは単一のベクトル PHI と、単一の観測値またはグループに対する行が含まれている XFUN を使用します。VFUN は、XFUN のすべての行に適用される単一行、または各行が XFUN の各行に対応する行列にすることができます。これは既定の設定です。 'SingleGroup' — modelfun は、データ内の 1 つのグループに対応する入力のみを受け入れるため、XFUN は呼び出しごとに 1 つのグループの X の行を含んでいなければなりません。モデルに応じて、PHI は、グループ全体に適用する 1 つの行、または観測ごとに 1 つの行をもつ行列です。VFUN は、1 つの行です。 'Full' — modelfun は、データ内の複数のパラメーター ベクトルと複数のグループの入力を受け入れます。PHI または VFUN は、XFUN のすべての行に適用する 1 つの行、または XFUN 内の行に対応する行をもつ行列です。このオプションを使用すると、modelfun への呼び出し回数が減ることでパフォーマンスが向上しますが、PHI または V で 1 次元拡張を実行するために modelfun が必要です。
CovParameterization	スケーリングされた共分散行列に対して内部的に使用されるパラメーター表現を指定します。選択肢は、コレスキー分解の場合は 'chol'、行列対数の場合は 'logm' です。既定の設定は 'logm' です。
CovPattern	r 行 r 列の logical 行列または数値行列 P を指定します。この行列は、変量効果の共分散行列 PSI のパターンを定義します。nlmefit は、PSI の対角に沿った分散と、P の非対角要素にある非ゼロによって指定される共分散を推定します。P の 0 の非対角要素に対応する共分散は 0 に制約されます。P がブロック対角行列の行と列の置換を指定しない場合は、nlmefit は、非ゼロ要素を必要に応じて P に追加します。P の既定値は、無相関の変量効果に対応する eye(r) です。 または、P は、変量効果のグループを指定する等しい値をもつ 1:r の値を含む 1 行 r 列のベクトルになります。この場合、nlmefit は、グループ内でのみ共分散を推定し、グループにまたがる共分散を 0 に制約します。
ParamTransform	P のパラメーター XB = ADESIGN * BETA + BDESIGN * B および PHI = f(XB) のそれぞれに対する変換関数 f() を指定する p 値のベクトル。ベクトルの各要素は、対応する PHI 値に変換を指定する次のいずれかの整数コードでなければなりません。 0: PHI = XB (すべてのパラメーターの既定値) 1: log(PHI) = XB 2: probit(PHI) = XB 3: logit(PHI) = XB
Options	statset によって返される形式の構造体。nlmefit は、次の statset パラメーターを使用します。 'DerivStep' — 有限差分の勾配の計算に使われる相対微分。スカラーまたは長さがモデル パラメーター p 数であるベクトル。既定の設定は eps^(1/3) です。 'Display' — 推定時の反復表示のレベル。選択肢は、以下のとおりです。 'off' (既定) — 何の情報も表示しません。 'final' — 最後の反復後に情報を表示します。 'iter' — 各反復で情報を表示します。 'FunValCheck' — NaN や Inf などの modelfun の無効な値をチェックします。選択肢は、'on' と 'off' です。既定の設定は 'on' です。 'MaxIter' — 許容される最大反復回数。既定の設定は 200 です。 'OutputFcn' — 関数ハンドルまたは空の配列 (既定値) をもつ cell 配列 @ を使用して指定された関数ハンドル。ソルバーは、各反復後にすべての出力関数を呼び出します。 'TolFun' — 対数尤度関数の終了許容誤差。既定の設定は 1e-4 です。 'TolX' — 推定された固定効果と変量効果の終了許容誤差。既定の設定は 1e-4 です。
OptimFun	尤度関数を最大化する推定プロセスの最適化関数。'fminsearch' (既定) または 'fminunc' として指定します。'fminunc' を使用するには、Optimization Toolbox™ が必要です。関数 fminsearch は、関数評価のみを使用する直接探索法を使用します。関数 fminunc (Optimization Toolbox) は勾配法を使用し、一般に、尤度関数を最大化する最適化問題ではこの関数の方が効率的です。

例

非線形混合効果モデル

5 本のオレンジのツリーの成長に関するデータを入力して表示します。

CIRC = [30 58 87 115 120 142 145;
        33 69 111 156 172 203 203;
        30 51 75 108 115 139 140;
        32 62 112 167 179 209 214;
        30 49 81 125 142 174 177];
time = [118 484 664 1004 1231 1372 1582];

h = plot(time,CIRC','o','LineWidth',2);
xlabel('Time (days)')
ylabel('Circumference (mm)')
title('{\bf Orange Tree Growth}')
legend([repmat('Tree ',5,1),num2str((1:5)')],...
       'Location','NW')
grid on
hold on

無名関数を使用して、ロジスティック成長モデルを指定します。

model = @(PHI,t)(PHI(:,1))./(1+exp(-(t-PHI(:,2))./PHI(:,3)));

既定の設定で、nlmefit を使用してモデルを当てはめます (つまり、各パラメーターが固定効果と変量効果の和であり、変量効果の間に相関がないと仮定します)。

TIME = repmat(time,5,1);
NUMS = repmat((1:5)',size(time));

beta0 = [100 100 100];
[beta1,PSI1,stats1] = nlmefit(TIME(:),CIRC(:),NUMS(:),...
                              [],model,beta0)

beta1 = 3×1

  191.3189
  723.7608
  346.2517

PSI1 = 3×3

  962.1534         0         0
         0    0.0000         0
         0         0  297.9882

stats1 = struct with fields:
           dfe: 28
          logl: -131.5457
           mse: 59.7882
          rmse: 7.9016
    errorparam: 7.7323
           aic: 277.0913
           bic: 274.3574
          covb: [3x3 double]
        sebeta: [15.2249 33.1579 26.8235]
          ires: [35x1 double]
          pres: [35x1 double]
         iwres: [35x1 double]
         pwres: [35x1 double]
         cwres: [35x1 double]

2 番目の変量効果の無視できる分散 PSI1(2,2) は、モデルを単純化するために取り除くことができることを示しています。

[beta2,PSI2,stats2,b2] = nlmefit(TIME(:),CIRC(:),...
    NUMS(:),[],model,beta0,'REParamsSelect',[1 3])

beta2 = 3×1

  191.3190
  723.7611
  346.2529

PSI2 = 2×2

  962.0448         0
         0  298.2191

stats2 = struct with fields:
           dfe: 29
          logl: -131.5457
           mse: 59.7879
          rmse: 7.7641
    errorparam: 7.7323
           aic: 275.0913
           bic: 272.7479
          covb: [3x3 double]
        sebeta: [15.2241 33.1578 26.8244]
          ires: [35x1 double]
          pres: [35x1 double]
         iwres: [35x1 double]
         pwres: [35x1 double]
         cwres: [35x1 double]

b2 = 2×5

  -28.5253   31.6061  -36.5070   39.0738   -5.6476
   10.0043   -0.7634    6.0086   -9.4638   -5.7858

対数尤度 logl は影響を受けず、赤池情報量基準とベイズ情報量基準 (aic と bic) はともに減少します。これはモデルから 2 番目の変量効果を棄却する決定をサポートしています。

beta2 の推定した固定効果と b2 のそれぞれのツリーに対して推定した変量効果を使用して、データからモデルをプロットします。

PHI = repmat(beta2,1,5) + ...          % Fixed effects
      [b2(1,:);zeros(1,5);b2(2,:)];    % Random effects

tplot = 0:0.1:1600;
for I = 1:5
  fitted_model=@(t)(PHI(1,I))./(1+exp(-(t-PHI(2,I))./ ... 
       PHI(3,I)));
  plot(tplot,fitted_model(tplot),'Color',h(I).Color, ...
	   'LineWidth',2)
end

参考文献

[1] Lindstrom, M. J., and D. M. Bates. “Nonlinear mixed-effects models for repeated measures data.” Biometrics. Vol. 46, 1990, pp. 673–687.

[2] Davidian, M., and D. M. Giltinan. Nonlinear Models for Repeated Measurements Data. New York: Chapman & Hall, 1995.

[3] Pinheiro, J. C., and D. M. Bates. “Approximations to the log-likelihood function in the nonlinear mixed-effects model.” Journal of Computational and Graphical Statistics. Vol. 4, 1995, pp. 12–35.

[4] Demidenko, E. Mixed Models: Theory and Applications. Hoboken, NJ: John Wiley & Sons, Inc., 2004.

バージョン履歴

R2008b で導入

参考

nlinfit | nlpredci | nlmefitsa

トピック

• 混合効果のモデル
• グループ化変数