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

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

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

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

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

Для примера мы можем сравнить подгонку со ссылкой probit к единице с ссылкой logit.

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

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