非線形ロジスティック回帰

問題の説明

ロジスティック回帰は特殊なタイプの回帰であり、何かの確率を他の変数の関数としてモデル化することが目標です。予測子ベクトルの集合 $$x_1,\dots ,x_N$$ について考えます。ここで、$$N$$ は観測値の個数、$$x_i$$ は $$i$$ 番目の観測値についての $$d$$ 個の予測子の値が含まれている列ベクトルです。$$x_i$$ の応答変数は $$Z_i$$ です。ここで $$Z_i$$ は、試行回数 $$n$$ および試行 $$i$$ の成功確率 $${\mu }_i$$ というパラメーターをもつ二項確率変数を表します。正規化された応答変数は $$Y_i=Z_i/n$$ であり、観測値 $$i$$ についての $$n$$ 回の試行における成功の比率です。応答 $$Y_i$$ は $$i=1,\dots ,N$$ について独立していると仮定します。各 $$i$$ について

$$Var(Y_i)=\frac{{\mu }_i(1-{\mu }_i)}{n}$$.

予測子変数 $$x_i$$ の関数として $${\mu }_i$$ をモデル化することを考えます。

線形ロジスティック回帰では、関数 fitglm を使用して次のように $${\mu }_i$$ を $$x_i$$ の関数としてモデル化できます。

$$\beta$$ は $$x_i$$ 内の予測子に乗算する係数の集合を表します。しかし、右辺で非線形関数が必要であると仮定します。

Statistics and Machine Learning Toolbox™ には、非線形回帰モデルを当てはめる関数はありますが、非線形ロジスティック回帰モデルを当てはめる関数はありません。この例では、ツールボックス関数を使用してこのようなモデルを当てはめる方法を示します。

直接最尤 (ML) 法

ML のアプローチでは、観測したデータの対数尤度を最大化します。尤度は、関数 binopdf で計算した二項確率 (または密度) 関数を使用して簡単に計算できます。

一般化最小二乗 (GLS)

関数 fitnlm を使用すると、非線形ロジスティック回帰モデルを推定できます。fitnlm は二項分布やリンク関数に対応していないので、はじめは意外に思うかもしれません。しかし、応答の平均と分散を指定した場合、fitnlm では一般化最小二乗 (GLS) をモデルの推定に使用できます。GLS が収束する場合、ML で解く場合と同じ一連の非線形方程式の解が求められ、$$\beta$$ が推定されます。GLS は、一般化線形モデルの疑似尤度推定にも使用できます。言い換えると、GLS と ML では同じまたは等価な解が得られるはずです。GLS による推定を実装するには、当てはめる非線形関数と二項分布の分散関数を用意します。

平均またはモデル関数

モデル関数は、$${\mu }_i$$ が $$\beta$$ によって変化する様子を説明します。fitnlm の場合、モデル関数は次のようになります。

重み関数

fitnlm は、名前と値のペアの引数 'Weights' を使用して関数ハンドルとして観測値の重みを受け入れます。このオプションを使用する場合、fitnlm では次のモデルを仮定します。

ここで、応答 $$Y_i$$ は独立していると仮定されています。$$w$$ は、$${\mu }_i$$ を受け入れて観測値の重みを返すカスタム関数ハンドルです。言い換えると、重みは応答の分散に反比例します。ロジスティック回帰モデルで使用する二項分布の場合、次のように重み関数を作成します。

fitnlm は、応答 $$Y_i$$ の分散を $${\sigma }^2/w({\mu }_i)$$ としてモデル化します。ここで $${\sigma }^2$$ は、GLS 推定には存在し、ロジスティック回帰モデルには存在しない追加パラメーターです。ただし、このパラメーターは通常は $$\beta$$ の推定に影響を与えず、$$Z_i$$ の値が二項分布に従うという仮定をチェックするための "分散パラメーター" を提供します。

直接 ML 法ではなく fitnlm を使用することのメリットとして、モデルの係数に対して仮説検定を実行し信頼区間を計算できます。

サンプル データの生成

ML と GLS の近似の違いを示すため、いくつかのサンプル データを生成します。$$x_i$$ は 1 次元であり、非線形ロジスティック回帰モデルにおける真の関数 $$f$$ は $$2\times 1$$ のベクトル $$\beta$$ によってパラメーター表現されたミカエリス・メンテン モデルであると仮定します。

myf = @(beta,x) beta(1)*x./(beta(2) + x);

$${\mu }_i$$ と $$\beta$$ の関係を指定するモデル関数を作成します。

mymodelfun = @(beta,x) 1./(1 + exp(-myf(beta,x)));

1 次元予測子のベクトルと真の係数のベクトル $$\beta$$ を作成します。

rng(300,'twister');
x    = linspace(-1,1,200)';
beta = [10;2];

各予測子について $${\mu }_i$$ 個の値が含まれているベクトルを計算します。

mu = mymodelfun(beta,x);

成功確率 $${\mu }_i$$ および試行回数 $$n$$ をもつ二項分布から応答 $$z_i$$ を生成します。

n = 50;
z = binornd(n,mu);

応答を正規化します。

y = z./n;

ML のアプローチ

ML のアプローチでは、負の対数尤度をベクトル $$\beta$$ の関数として定義してから、fminsearch などの最適化関数で最小化します。$$\beta$$ の開始値として beta0 を指定します。

mynegloglik = @(beta) -sum(log(binopdf(z,n,mymodelfun(beta,x))));
beta0 = [3;3];
opts = optimset('fminsearch');
opts.MaxFunEvals = Inf;
opts.MaxIter = 10000;
betaHatML = fminsearch(mynegloglik,beta0,opts)

betaHatML = 2×1

    9.9259
    1.9720

betaHatML の推定係数は、真の値である [10;2] に近くなっています。

GLS のアプローチ

GLS のアプローチでは、前述した fitnlm の重み関数を作成します。

wfun = @(xx) n./(xx.*(1-xx));

平均と重みのカスタム関数を指定して fitnlm を呼び出します。$$\beta$$ の開始値として beta0 を指定します。

nlm = fitnlm(x,y,mymodelfun,beta0,'Weights',wfun)

nlm = 
Nonlinear regression model:
    y ~ F(beta,x)

Estimated Coefficients:
          Estimate      SE       tStat       pValue  
          ________    _______    ______    __________

    b1     9.926      0.83135     11.94     4.193e-25
    b2     1.972      0.16438    11.996    2.8182e-25


Number of observations: 200, Error degrees of freedom: 198
Root Mean Squared Error: 1.16
R-Squared: 0.995,  Adjusted R-Squared 0.995
F-statistic vs. zero model: 1.88e+04, p-value = 2.04e-226

当てはめた NonLinearModel オブジェクト nlm から $$\beta$$ の推定を取得します。

betaHatGLS = nlm.Coefficients.Estimate

betaHatGLS = 2×1

    9.9260
    1.9720

ML 法の場合と同じように、betaHatGLS の推定係数は真の値である [10;2] に近くなっています。$${\beta }_1$$ と $${\beta }_2$$ の両方について p 値が小さいので、どちらの係数も $$0$$ とは有意に異なることがわかります。

ML のアプローチと GLS のアプローチの比較

$$\beta$$ の推定を比較します。

max(abs(betaHatML - betaHatGLS))

ans = 1.1460e-05

ML を使用した近似値と GLS を使用した近似値を比較します。

yHatML  = mymodelfun(betaHatML ,x);
yHatGLS = mymodelfun(betaHatGLS,x);
max(abs(yHatML - yHatGLS))

ans = 1.2746e-07

ML のアプローチと GLS のアプローチでは、同様の解が得られます。

ML と GLS を使用した近似値のプロット

figure
plot(x,y,'g','LineWidth',1)
hold on
plot(x,yHatML ,'b'  ,'LineWidth',1)
plot(x,yHatGLS,'m--','LineWidth',1)
legend('Data','ML','GLS','Location','Best')
xlabel('x')
ylabel('y and fitted values')
title('Data y along with ML and GLS fits.')

ML と GLS では、同様の近似値が得られます。

ML と GLS を使用して推定した非線形関数のプロット

$$f(x_i,\beta )$$ について真のモデルをプロットします。$$\beta ={\beta }_0$$ を使用した $$f(x_i,\beta )$$ の初期推定についてのプロットと、ML ベースおよび GLS ベースの $$f(x_i,\beta )$$ の推定のプロットを追加します。

figure
plot(x,myf(beta,x),'r','LineWidth',1)
hold on
plot(x,myf(beta0,x),'k','LineWidth',1)
plot(x,myf(betaHatML,x),'c--','LineWidth',1)
plot(x,myf(betaHatGLS,x),'b-.','LineWidth',1)
legend('True f','Initial f','Estimated f with ML', ...
    'Estimated f with GLS','Location','Best')
xlabel('x')
ylabel('True and estimated f')
title('Comparison of true f with estimated f using ML and GLS.')

ML 法と GLS 法の両方を使用して推定された非線形関数 $$f$$ は、真の非線形関数 $$f$$ に近くなっています。同様の手法を使用して、非線形ポアソン回帰など他の非線形一般化線形モデルを当てはめることができます。