ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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 は、カテゴリカル変数、数値ベクトル、グループ名の行をもつ文字行列、または文字列のセル配列になります。グループ化変数の詳細は、「グループ化変数」を参照してください。

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

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

yfit = modelfun(PHI,XFUN,VFUN)

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

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

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

    • k = 1 (XFUNX の単一行の場合)

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

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

  • 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 — 誤差分散モデルの推定パラメーター

  • aicaic = -2 * logl + 2 * numParam の式で計算される赤池情報量基準。ここで、numParam は変量効果の共分散行列の自由度を含めた近似パラメーターの数、固定効果の数、および誤差モデルのパラメーターの数、loglstats 構造体のフィールドです。

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

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

    • numParamlogl は、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 を指定します。

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

パラメーター
FEParamsSelect

1p のインデックスの数値ベクトルまたは 1 行 p 列の論理ベクトルとして与えられた、固定効果を含むパラメーターのベクトル PHI の要素を指定するベクトル。q が指定した要素数である場合、モデルは q 個の固定効果を含みます。

FEConstDesign

p 行 q 列の計画行列 ADESIGN。ここで、ADESIGN*betaPHI の p 要素の固定要素です。

FEGroupDesign

m 個のグループごとに異なる p 行 q 列の固定効果の計画行列を指定する p x q x m 配列。

FEObsDesign

n 個の観測ごとに異なる p 行 q 列の固定効果の計画行列を指定する p x q x n 配列。

REParamsSelect

1p のインデックスの数値ベクトルまたは 1 行 p 列の論理ベクトルとして与えられた、変量効果を含むパラメーターのベクトル PHI の要素を指定するベクトル。モデルは r 個の変量効果を含みます。ここで、r は、指定した要素数です。

REConstDesign

p 行 r 列の計画行列 BDESIGN。ここで、BDESIGN*BPHI の p 要素のランダム要素です。

REGroupDesign

m 個のグループごとに異なる p 行 r 列の変量効果の計画行列を指定する p x r x m 配列。

REObsDesign

n 個の観測ごとに異なる p 行 r 列の変量効果の計画行列を指定する p x r x n 配列。

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

Parameter
RefineBeta0

nlmefit が変量効果のない modelfun を近似した後、beta0beta で置き換えることで、beta0 の初期調整を行うかどうかを指定します。選択肢は、'on''off' です。既定値は 'on' です。

ErrorModel

誤差項の形式を指定する文字列。既定の設定は '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'betaB の現在の条件推定で線形混合効果モデルに対する尤度を使用します。これは既定値です。

  • 'RELME'betaB の現在の条件推定で線形混合効果モデルに対する制限付き尤度を使用します。

  • 'FO' — 変量効果のない 1 次ラブラシアン近似。

  • 'FOCE'B の条件推定での 1 次ラブラシアン近似。

Vectorization

modelfun への PHIXFUNVFUN の入力引数の許容可能なサイズを指定します。選択肢は、以下のとおりです。

  • 'SinglePhi'modelfun は一度に 1 つのモデル パラメーター セットのみを受け入れるため、PHI は呼び出しごとに 1 つの行ベクトルでなければなりません。nlmefit は、必要に応じて、1 つの PHI ベクトルと、一度に 1 つの観測またはグループに対する行を含む XFUN を使用して modelfun をループ内で呼び出します。VFUN は、XFUN のすべての行に適用する 1 つの行、または 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 列の論理行列または数値行列 P を指定します。この行列は、変量効果の共分散行列 PSI のパターンを定義します。nlmefit は、PSI の対角に沿った分散と、P の非対角要素にある非ゼロによって指定される共分散を推定します。P の 0 の非対角要素に対応する共分散は 0 に制約されます。P がブロック対角行列の行と列の置換を指定しない場合は、nlmefit は、非ゼロ要素を必要に応じて P に追加します。P の既定値は、無相関の変量効果に対応する eye(r) です。

または、P は、変量効果のグループを指定する等しい値をもつ 1:r の値を含む 1 行 r 列のベクトルになります。この場合、nlmefit は、グループ内でのみ共分散を推定し、グループにまたがる共分散を 0 に制約します。

ParamTransform

P パラメーターごとに変換関数 f() を指定する p 値のベクトル XB = ADESIGN*BETA + BDESIGN*B PHI = f(XB)。ベクトルの各要素は、対応する 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'NaNInf などの modelfun の無効な値をチェックします。選択肢は、'on''off' です。既定の設定は 'on' です。

  • 'MaxIter' — 許容される最大反復回数。既定の設定は 200 です。

  • 'OutputFcn' — 関数ハンドルまたは空の配列 (既定値) をもつセル配列 @ を使用して指定された関数ハンドル。ソルバーは、各反復後にすべての出力関数を呼び出します。

  • 'TolFun' — 対数尤度関数の終了許容誤差。既定の設定は 1e-4 です。

  • 'TolX' — 推定された固定効果と変量効果の終了許容誤差。既定値は 1e-4 です。

OptimFun

尤度の最大化で使用する最適化関数を指定します。選択肢は、fminsearch を使用するには 'fminsearch'fminunc を使用するには 'fminunc' です。既定の設定は 'fminsearch' です。Optimization Toolbox™ ソフトウェアがインストールされている場合は、'fminunc' のみ指定できます。

すべて展開する

非線形混合効果モデル

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 =

  191.3189
  723.7608
  346.2517


PSI1 =

  962.1533         0         0
         0    0.0000         0
         0         0  297.9880


stats1 = 

           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 =

  191.3187
  723.7603
  346.2483


PSI2 =

  962.5842         0
         0  297.2055


stats2 = 

           dfe: 29
          logl: -131.5457
           mse: 59.7884
          rmse: 7.7647
    errorparam: 7.7323
           aic: 275.0913
           bic: 272.7480
          covb: [3x3 double]
        sebeta: [15.2277 33.1583 26.8205]
          ires: [35x1 double]
          pres: [35x1 double]
         iwres: [35x1 double]
         pwres: [35x1 double]
         cwres: [35x1 double]


b2 =

  -28.5282   31.6066  -36.5090   39.0769   -5.6463
    9.9752   -0.7600    5.9907   -9.4372   -5.7687

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

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

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

colors = get(h,'Color');
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',colors{I}, ...
			'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.

この情報は役に立ちましたか?