Main Content

일반화 선형 모델을 사용하여 데이터 피팅하기

이 예제에서는 glmfitglmval을 사용하여 일반화 선형 모델을 피팅하고 평가하는 방법을 보여줍니다. 보통 선형 회귀를 사용하여 직선, 즉 해당 모수에서 선형인 함수를 정규분포된 오차를 갖는 데이터에 피팅할 수 있습니다. 이것이 가장 일반적으로 사용되는 회귀 모델입니다. 그러나, 항상 현실적인 모델은 아닙니다. 일반화 선형 모델은 선형 모델을 두 가지 관점에서 확장합니다. 첫째, 연결 함수 추가를 통해 모수의 선형성에 대한 가정이 완화됩니다. 둘째, 정규분포와 다른 오차 분포를 모델링할 수 있습니다.

일반화 선형 모델

회귀 모델은 하나 이상의 예측 변수(대개 x1, x2 등으로 표시됨)에 대해 응답 변수(일반적으로 대개 y로 표시됨)의 분포를 정의합니다. 가장 일반적으로 사용되는 회귀 모델인 보통 선형 회귀는 y를 평균이 예측 변수의 선형 함수(b0 + b1*x1 + ...)이고 분산이 상수인 정규 확률 변수로 모델링합니다. 단일 예측 변수 x의 가장 단순한 사례에서는 모델을 각 점에 대한 가우스 분포를 갖는 직선으로 나타낼 수 있습니다.

mu = @(x) -1.9+.23*x;
x = 5:.1:15;
yhat = mu(x);
dy = -3.5:.1:3.5; sz = size(dy); k = (length(dy)+1)/2;
x1 =  7*ones(sz); y1 = mu(x1)+dy; z1 = normpdf(y1,mu(x1),1);
x2 = 10*ones(sz); y2 = mu(x2)+dy; z2 = normpdf(y2,mu(x2),1);
x3 = 13*ones(sz); y3 = mu(x3)+dy; z3 = normpdf(y3,mu(x3),1);
plot3(x,yhat,zeros(size(x)),'b-', ...
      x1,y1,z1,'r-', x1([k k]),y1([k k]),[0 z1(k)],'r:', ...
      x2,y2,z2,'r-', x2([k k]),y2([k k]),[0 z2(k)],'r:', ...
      x3,y3,z3,'r-', x3([k k]),y3([k k]),[0 z3(k)],'r:');
zlim([0 1]);
xlabel('X'); ylabel('Y'); zlabel('Probability density');
grid on; view([-45 45]);

일반화 선형 모델에서 응답 변수의 평균은 예측 변수의 선형 함수에 대한 단조 비선형 변환(g(b0 + b1*x1 + ...))으로 모델링됩니다. 변환 g의 역함수를 "연결" 함수라고 합니다. 로짓(시그모이드) 연결과 로그 연결이 그러한 예입니다. 또한, y는 이항분포 또는 푸아송 분포와 같은 비정규분포를 가질 수 있습니다. 예를 들어, 로그 링크와 단일 예측 변수 x를 갖는 푸아송 회귀는 각 점에 대한 푸아송 분포를 갖는 지수 곡선으로 나타낼 수 있습니다.

mu = @(x) exp(-1.9+.23*x);
x = 5:.1:15;
yhat = mu(x);
x1 =  7*ones(1,5);  y1 = 0:4; z1 = poisspdf(y1,mu(x1));
x2 = 10*ones(1,7); y2 = 0:6; z2 = poisspdf(y2,mu(x2));
x3 = 13*ones(1,9); y3 = 0:8; z3 = poisspdf(y3,mu(x3));
plot3(x,yhat,zeros(size(x)),'b-', ...
      [x1; x1],[y1; y1],[z1; zeros(size(y1))],'r-', x1,y1,z1,'r.', ...
      [x2; x2],[y2; y2],[z2; zeros(size(y2))],'r-', x2,y2,z2,'r.', ...
      [x3; x3],[y3; y3],[z3; zeros(size(y3))],'r-', x3,y3,z3,'r.');
zlim([0 1]);
xlabel('X'); ylabel('Y'); zlabel('Probability');
grid on; view([-45 45]);

로지스틱 회귀 피팅하기

이 예제에서는 다양한 중량의 차량에서 주행 거리 검정에 실패한 차량의 비율을 모델링하는 데 도움이 되는 실험을 다룹니다. 데이터는 중량, 검정된 차량 수, 실패한 차량 수에 대한 관측값을 포함합니다.

% A set of car weights
weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
% The number of cars tested at each weight
tested = [48 42 31 34 31 21 23 23 21 16 17 21]';
% The number of cars failing the test at each weight
failed = [1 2 0 3 8 8 14 17 19 15 17 21]';
% The proportion of cars failing for each weight
proportion = failed ./ tested;

plot(weight,proportion,'s')
xlabel('Weight'); ylabel('Proportion');

이 그래프는 중량에 대한 함수로 실패하는 차량의 비율을 보여주는 플롯입니다. 중량과 함께 증가하는 확률 모수 P를 사용하여 실패 수가 이항분포에서 왔다고 가정하는 것이 합리적입니다. 그러나 P가 얼마나 정확히 중량에 종속되어야 할까요?

직선을 이러한 데이터에 피팅해 볼 수 있습니다.

linearCoef = polyfit(weight,proportion,1);
linearFit = polyval(linearCoef,weight);
plot(weight,proportion,'s', weight,linearFit,'r-', [2000 4500],[0 0],'k:', [2000 4500],[1 1],'k:')
xlabel('Weight'); ylabel('Proportion');

이 선형 피팅에는 두 가지 문제가 있습니다.

1) 선은 0보다 작고 1보다 큰 비율을 예측합니다.

2) 비율이 유계여야 하므로 정규분포되지 않습니다. 이는 단순한 선형 회귀 모델을 피팅하는 데 필요한 가정 중 하나를 위반합니다.

고차 다항식을 사용하는 것이 도움이 되는 것처럼 보일 수 있습니다.

[cubicCoef,stats,ctr] = polyfit(weight,proportion,3);
cubicFit = polyval(cubicCoef,weight,[],ctr);
plot(weight,proportion,'s', weight,cubicFit,'r-', [2000 4500],[0 0],'k:', [2000 4500],[1 1],'k:')
xlabel('Weight'); ylabel('Proportion');

그러나, 이 피팅에도 여전히 유사한 문제가 있습니다. 이 그래프는 중량이 4000을 넘어서면서 피팅된 비율이 감소함을 보여줍니다. 실제로, 이 비율은 중량 값이 더 커지면 음수가 됩니다. 당연히, 정규분포 가정을 계속 위반합니다.

대신, glmfit을 사용하여 로지스틱 회귀 모델을 피팅하는 것이 더 나은 방법입니다. 로지스틱 회귀는 일반화 선형 모델의 특수한 사례이며 다음 두 가지 이유로 이 데이터에 대해서는 선형 회귀보다 더 적합합니다. 첫째, 이항분포에 적합한 피팅 방법을 사용합니다. 둘째, 로지스틱 연결이 예측된 비율을 범위 [0,1]로 제한합니다.

로지스틱 회귀를 위해 예측 변수 행렬과 함께 하나의 행렬을 더 지정합니다. 이 행렬의 한 열은 실패 수를 포함하고 다른 열은 검정된 수를 포함합니다. 또한, 이항분포와 로짓 연결도 지정합니다.

[logitCoef,dev] = glmfit(weight,[failed tested],'binomial','logit');
logitFit = glmval(logitCoef,weight,'logit');
plot(weight,proportion,'bs', weight,logitFit,'r-');
xlabel('Weight'); ylabel('Proportion');

이 플롯에 나타난 것처럼 피팅된 비율은 중량이 작아지면 0에 가까워지고 중량이 커지면 1에 가까워집니다.

모델 진단

glmfit 함수는 여러 출력값을 제공합니다. 이를 통해 피팅을 검토하고 모델을 검정할 수 있습니다. 예를 들어, 두 모델의 이탈도 값을 비교하여 제곱 항이 피팅을 크게 향상시키는지를 확인할 수 있습니다.

[logitCoef2,dev2] = glmfit([weight weight.^2],[failed tested],'binomial','logit');
pval = 1 - chi2cdf(dev-dev2,1)
pval =

    0.4019

p-값이 클 경우 이 데이터에 대한 2차 항이 피팅을 크게 향상시키지 않는다는 것을 나타냅니다. 두 피팅에 대한 플롯은 두 피팅에 차이가 거의 없음을 보여줍니다.

logitFit2 = glmval(logitCoef2,[weight weight.^2],'logit');
plot(weight,proportion,'bs', weight,logitFit,'r-', weight,logitFit2,'g-');
legend('Data','Linear Terms','Linear and Quadratic Terms','Location','northwest');

피팅의 적합도를 확인하기 위해 피어슨 잔차에 대한 확률 플롯을 살펴볼 수도 있습니다. 이 잔차는 모델이 데이터에 적절한 피팅인 경우 표준 정규분포에 거의 가까운 분포를 가지도록 정규분포됩니다. (이 표준화가 적용되지 않을 경우 잔차가 다른 분산을 가질 수 있습니다.)

[logitCoef,dev,stats] = glmfit(weight,[failed tested],'binomial','logit');
normplot(stats.residp);

이 잔차 플롯은 정규분포를 훌륭히 따름을 보여줍니다.

모델 예측 실행하기

모델에 만족하면, 이 모델을 사용하여 신뢰한계 계산을 포함하여 예측을 수행할 수 있습니다. 여기서는 검정되는 차량 100대 중에 4개 중량 각각에서 주행 거리 검정에 실패할 것으로 예상되는 수를 예측하겠습니다.

weightPred = 2500:500:4000;
[failedPred,dlo,dhi] = glmval(logitCoef,weightPred,'logit',stats,.95,100);
errorbar(weightPred,failedPred,dlo,dhi,':');

이항 모델에 대한 연결 함수

glmfit이 지원하는 5개 분포 각각에는 정준(디폴트) 연결 함수가 있습니다. 이항분포의 경우 정준 연결은 로짓입니다. 그러나, 이항 모델에 적합한 다른 연결 함수가 3개 더 있습니다. 이 4개 연결 함수 모두 평균 응답 변수를 구간 [0, 1] 내에 유지합니다.

eta = -5:.1:5;
plot(eta,1 ./ (1 + exp(-eta)),'-', eta,normcdf(eta), '-', ...
     eta,1 - exp(-exp(eta)),'-', eta,exp(-exp(eta)),'-');
xlabel('Linear function of predictors'); ylabel('Predicted mean response');
legend('logit','probit','complementary log-log','log-log','location','east');

예를 들어, 프로빗 연결 함수를 사용한 피팅과 로짓 연결 함수를 사용한 피팅을 비교할 수 있습니다.

probitCoef = glmfit(weight,[failed tested],'binomial','probit');
probitFit = glmval(probitCoef,weight,'probit');
plot(weight,proportion,'bs', weight,logitFit,'r-', weight,probitFit,'g-');
legend('Data','Logit model','Probit model','Location','northwest');

데이터가 이 4개 연결 함수 간을 구분하는 것은 대개 어려우며 이론적 근거에 기반하여 선택하는 경우가 많습니다.