Main Content

가중 비선형 회귀

이 예제에서는 상수가 아닌 오차 분산을 갖는 데이터에 비선형 회귀 모델을 피팅하는 방법을 보여줍니다.

정규 비선형 최소제곱 알고리즘은 계측 오차 모두 동일한 분산을 갖는 경우에 적합합니다. 이러한 가정이 충족되지 않으면 가중치 피팅을 사용하는 것이 적합합니다. 이 예제에서는 fitnlm 함수에 가중치를 사용하는 방법을 보여줍니다.

피팅을 위한 데이터 및 모델

산업 및 가정용 폐기물로 인한 수질 오염을 연구하기 위해 수집된 데이터를 사용하겠습니다. 이 데이터는 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')

Figure contains an axes object. The axes object contains an object of type line.

처음 두 관측값이 나머지 관측값보다 낮은 정밀도로 측정되었음을 알고 있다고 가정하겠습니다. 예를 들어, 처음 두 관측값이 다른 계측기로 측정되었을 수 있습니다. 데이터에 가중치를 적용하는 또 다른 일반적인 경우는 기록된 각 관측값이 실제로 동일한 x 값에서 측정한 여러 값의 평균일 때입니다. 여기 데이터에서는 처음 두 값이 단일 원시 측정값을 나타내는 반면, 나머지 4개 값은 각각 5개 원시 측정값의 평균에 해당한다고 가정하겠습니다. 그러면 각 관측값에 사용된 측정값 개수를 기준으로 가중치를 적용하는 것이 적합할 수 있습니다.

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

가중치 없이 모델 피팅하기

이 데이터에 피팅하는 모델은 x가 커지면 평평해지는 스케일링된 지수 곡선입니다.

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

거친 시각적 피팅을 바탕으로 점을 통과하여 그려진 곡선이 x = 15 근방의 약 240인 값에서 평평해지는 것처럼 보입니다. 따라서 240을 b1의 시작 값으로 사용하고, e^(-.5*15)이 1보다 작으므로 0.5를 b2의 시작 값으로 사용하겠습니다.

start = [240; .5];

계측 오차를 무시하는 경우 위험성은 피팅이 부정확한 측정값으로 인해 과도하게 영향을 받을 수 있으므로 정확하게 알려진 측정값에 대한 양호한 피팅을 제공하지 않을 수 있다는 점입니다. 가중치 없이 데이터를 피팅하고 이를 점과 비교해 보겠습니다.

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

Figure contains an axes object. The axes object contains 2 objects of type line.

피팅된 곡선이 처음 두 점에는 가깝게 나왔지만, 다른 점의 추세는 놓치는 것처럼 보입니다.

가중치를 사용하여 모델 피팅하기

가중치를 사용하여 피팅을 반복해 보겠습니다.

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')

Figure contains an axes object. The axes object contains 3 objects of type line.

이 사례의 추정된 모집단 표준편차는 가중치, 즉 측정 정밀도 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')

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent Data, Weighted fit, 95% Confidence Limits.

비중을 낮춘 두 개의 점이 나머지 점만큼 곡선에 잘 피팅되지 않았습니다. 가중치 피팅에 대해 예상한 결과 그대로입니다.

또한, 지정된 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')

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent Data, Weighted fit, 95% Prediction Limits.

가중치의 절대 스케일은 실제로 모수 추정값에 영향을 미치지 않습니다. 어떠한 상수로 가중치를 다시 스케일링하든 추정값은 동일합니다. 그러나, 신뢰한계가 가중치 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')

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent Data, Weighted fit, 95% Prediction Limits.

predict 함수는 MSE*(1/W(i))에 의해 관측값 i에서의 오차 분산을 추정합니다. 여기서 MSE는 평균 제곱 오차입니다. 그러므로 신뢰구간이 좁아집니다.

잔차 분석

데이터와 피팅을 플로팅하는 것 외에 예측 변수에 대해 피팅의 잔차를 플로팅하여 모델 관련 문제를 진단합니다. 잔차는 독립적이며 똑같이 분산되지만(i.i.d.) 분산은 가중치의 역에 비례한 것처럼 보여야 합니다. 표준화 잔차를 플로팅하여 값이 i.i.d이고 분산이 동일한지 확인합니다. 표준화 잔차는 잔차의 추정된 표준편차로 나눈 원시 잔차입니다.

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

Figure contains an axes object. The axes object contains an object of type line.

이 잔차 플롯에는 규칙적인 패턴을 보여주는 증거가 있습니다. 마지막 4개 잔차가 어떻게 선형 추세를 가지는지 살펴보십시오. 그러면 모델이 x가 증가함에 따라 충분히 빠르게 증가하지 않을 수 있다는 것을 알 수 있습니다. 또한, 잔차의 크기가 x가 증가할수록 감소하는 경향이 있으며, 이는 계측 오차가 x에 종속될 수 있음을 나타냅니다. 이는 조사할 가치가 있지만, 데이터 점이 극히 적기 때문에 이러한 분명한 패턴에 유의성을 부여하기가 어렵습니다.

참고 항목

| | |