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

Модель регрессии задает распределение переменной отклика (часто в общем обозначенный, когда 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');

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

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

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

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