Обобщенные линейные модели

Регрессионная модель определяет распределение переменной отклика (часто обобщенно обозначаемой как y) в терминах одной или нескольких переменных-предикторов (часто обозначаемых x1, x2 и т. д.). Наиболее часто используемая регрессионная модель, обычная линейная регрессия, моделирует 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 известно как функция «связи». Примеры включают ссылку logit (sigmoid) и ссылку log. Кроме того, 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');

Как показывает этот график, аппроксимированные пропорции асимптантны нулю и единице по мере того, как вес становится малым или большим.

Диагностика модели

glmfit функция предоставляет ряд выходных данных для проверки соответствия и тестирования модели. Например, мы можем сравнить значения отклонения для двух моделей, чтобы определить, будет ли квадрат значительно улучшить вписывание.

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

pval =

    0.4019

Большое значение p указывает на то, что для этих данных квадратичный член существенно не улучшает подгонку. График двух посадок показывает, что в посадках мало различий.

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 протестированных, которые провалили бы испытание пробега на каждом из четырех весов.

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

Функции связи для биномиальных моделей

Для каждого из пяти распределений, которые glmfit поддерживает, существует каноническая (по умолчанию) функция связи. Для биномиального распределения каноническим звеном является логит. Тем не менее, есть также три других ссылки, которые разумны для биномиальных моделей. Все четыре поддерживают средний отклик в интервале [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');

Данные часто трудно различить между этими четырьмя функциями связи, и выбор часто делается на теоретических основаниях.