重み付き非線形回帰

近似のためのデータとモデル

収集したデータを使用して、産業廃棄物と家庭廃棄物に起因する水質汚染について研究します。これらのデータは、Box, G.P., W.G. Hunter, and J.S. Hunter, Statistics for Experimenters (Wiley, 1978, pp. 483-487) で詳しく説明されています。応答変数は、生化学的酸素要求量 (mg/l) で、予想変数は培養時間 (日数) です。

x = [1 2 3 5 7 10]';
y = [109 149 149 191 213 224]';

plot(x,y,'ko')
xlabel('Incubation (days), x') 
ylabel('Biochemical oxygen demand (mg/l), y')

最初 2 回の観測は、残りの観測よりも低い精度で行われたことがわかっていると仮定します。たとえば、異なる器具を使用して行った場合などです。データに重みを付けるもう 1 つの理由は、記録された各観測値が実際には同じ値 x を使用して行った複数の測定の平均であるためです。ここで使用するデータでは、最初の 2 つの値が 1 回の生の測定値を表し、残りの 4 つの値がそれぞれ 5 回の生の測定値の平均であるとします。すると、各観測で使用した測定回数によって重みを付けることが適切です。

w = [1 1 5 5 5 5]';

重みを使用しないモデルの当てはめ

これらのデータに当てはめるモデルは、x が大きくなると共に平坦になるスケーリングされた指数曲線です。

modelFun = @(b,x) b(1).*(1-exp(-b(2).*x));

おおよその視覚的な近似に基づくと、点を経由して描画された曲線は 240 付近の値 (x = 15 付近) で横ばい状態になります。したがって、240 を b1 の開始値として使用し、e^(-.5*15) は 1 に比べて小さいので、.5 を b2 の開始値として使用します。

start = [240; .5];

測定誤差を無視すると、近似が不正確な測定の影響を過度に受けてしまう危険性があり、さらにこれが原因で既知の正確な測定への近似ができなくなる可能性があります。重みを使用せずにデータを近似し、各点と比較してみましょう。

nlm = fitnlm(x,y,modelFun,start);
xx = linspace(0,12)';
line(xx,predict(nlm,xx),'linestyle','--','color','k')

近似曲線は最初の 2 点には引き寄せられていますが、他の点のトレンドには従っていないようです。

重みを使用したモデルの当てはめ

重みを使用した近似を繰り返してみましょう。

wnlm = fitnlm(x,y,modelFun,start,'Weight',w)

wnlm = 
Nonlinear regression model:
    y ~ b1*(1 - exp( - b2*x))

Estimated Coefficients:
          Estimate       SE       tStat       pValue  
          ________    ________    ______    __________

    b1     225.17         10.7    21.045    3.0134e-05
    b2    0.40078     0.064296    6.2333     0.0033745


Number of observations: 6, Error degrees of freedom: 4
Root Mean Squared Error: 24
R-Squared: 0.908,  Adjusted R-Squared 0.885
F-statistic vs. zero model: 696, p-value = 8.2e-06

line(xx,predict(wnlm,xx),'color','b')

この場合、推定母標準偏差は、重みまたは測定精度 1 の "標準" 観測値の平均誤差を示します。

wnlm.RMSE

ans = 
24.0096

どのような分析でも、モデル近似の精度の推定は重要な部分です。係数表示にはパラメーターの標準誤差が表示されますが、信頼区間を計算することもできます。

coefCI(wnlm)

ans = 2×2

  195.4650  254.8788
    0.2223    0.5793

応答曲線の推定

次に、近似した応答値と、それに対する信頼区間を計算します。既定では、この幅は予測値の点単位の信頼限界に対するものですが、ここでは曲線全体の同時区間を計算します。

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true);
plot(x,y,'ko',xx,ypred,'b-',xx,ypredci,'r:')
xlabel('x') 
ylabel('y')
ylim([-150 350])
legend({'Data','Weighted fit','95% Confidence Limits'}, ...
    'location','SouthEast')

曲部では、重みを減らした 2 つの点は、残りの点ほど正確に近似されないことに注意してください。つまり、重み付き近似には期待されるほどの効果がないということです。

指定値 x で将来の観測の予測区間を推定することもできます。これらの区間では、重み (測定精度) が 1 と実際に仮定されます。

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true, ...
    'Prediction','observation');
plot(x,y,'ko',xx,ypred,'b-',xx,ypredci,'r:')
xlabel('x') 
ylabel('y')
ylim([-150 350])
legend({'Data','Weighted fit','95% Prediction Limits'}, ...
    'location','SouthEast')

重みの絶対スケールは、パラメーターの推定には影響しません。重みをどの定数で再スケーリングしても、推定結果は同じです。ただし、信頼限界には影響を与えます。これは、重み 1 の観測値を信頼限界が表しているためです。ここでは、信頼限界と比較すると重みが大きい点は、近似直線に近すぎることがわかります。

ここでは、このプロットの最後の 4 点と同様に、5 つの測定値の平均に基づく新しい観測を行います。関数 predict の名前と値の引数 Weights を使用して観測値の重みを指定します。

[new_ypred,new_ypredci] = predict(wnlm,xx,'Simultaneous',true, ...
    'Prediction','observation','Weights',5*ones(size(xx)));
plot(x,y,'ko',xx,new_ypred,'b-',xx,new_ypredci,'r:')
xlabel('x') 
ylabel('y')
ylim([-150 350])
legend({'Data','Weighted fit','95% Prediction Limits'}, ...
    'location','SouthEast')

関数 predict は、観測値 i における誤差分散を MSE*(1/W(i)) により推定します。MSE は平均二乗誤差です。したがって、信頼区間が狭くなります。

残差分析

データと近似のプロットに加え、予測子に対する近似からの残差もプロットし、モデルの問題を診断します。残差は独立同一分布で表示されますが、重みの逆数に比例する分散があります。標準化された残差をプロットして、値が同じ分散をもつ独立同一分布になることを確認します。標準化された残差は、生の残差を推定標準偏差で除算した値です。

r = wnlm.Residuals.Standardized;
plot(x,r,'b^')
xlabel('x')
ylabel('Standardized Residuals')

この残差プロットには、系統的なパターンが見られます。最後の 4 つの残差に線形のトレンドがあることに注意してください。この線形のトレンドとは、モデルが x の増大するスピードほど速く増大しないことを示しています。さらに、x が増加すると残差の大きさが減少することは、測定誤差が x に依存する可能性を示しています。これを調査する価値はありますが、データ点が少なすぎるため、この明らかなパターンに意味を持たせることは困難です。