В этом примере показано, как применять частичную регрессию наименьших квадратов (PLSR) и регрессию основных компонентов (PCR), и обсуждается эффективность этих двух методов. PLSR и PCR являются методами моделирования переменной ответа, когда существует большое количество переменных предиктора, и эти предикторы сильно коррелируют или даже коллинеарны. Оба метода строят новые переменные предиктора, известные как компоненты, как линейные комбинации исходных переменных предиктора, но они строят эти компоненты по-разному. ПЦР создает компоненты для объяснения наблюдаемой изменчивости переменных предиктора, вообще не рассматривая переменную ответа. С другой стороны, PLSR действительно учитывает переменную отклика и, следовательно, часто приводит к моделям, которые способны совместить переменную отклика с меньшим количеством компонентов. Будет ли это в конечном итоге преобразовано в более благоразумную модель с точки зрения ее практического использования, зависит от контекста.
Загрузить набор данных, содержащий спектральные интенсивности 60 образцов бензина на 401 длине волны и их октановые числа. Эти данные описаны в работе Kalivas, John H., «Два набора данных спектров ближнего инфракрасного излучения», Chemometrics and Intelligent Laboratory Systems, v.37 (1997) pp.255-259.
load spectra whos NIR octane
Name Size Bytes Class Attributes NIR 60x401 192480 double octane 60x1 480 double
[dummy,h] = sort(octane); oldorder = get(gcf,'DefaultAxesColorOrder'); set(gcf,'DefaultAxesColorOrder',jet(60)); plot3(repmat(1:401,60,1)',repmat(octane(h),1,401)',NIR(h,:)'); set(gcf,'DefaultAxesColorOrder',oldorder); xlabel('Wavelength Index'); ylabel('Octane'); axis('tight'); grid on
Используйте plsregress для размещения модели PLSR с десятью компонентами PLS и одним откликом.
X = NIR;
y = octane;
[n,p] = size(X);
[Xloadings,Yloadings,Xscores,Yscores,betaPLS10,PLSPctVar] = plsregress(...
X,y,10);
Десять компонентов могут быть больше, чем потребуется для адекватного подбора данных, но диагностика из этого подбора может быть использована для выбора более простой модели с меньшим количеством компонентов. Например, одним из быстрых способов выбора количества компонентов является построение графика процента дисперсии, объясненного в переменной отклика, как функции количества компонентов.
plot(1:10,cumsum(100*PLSPctVar(2,:)),'-bo'); xlabel('Number of PLS components'); ylabel('Percent Variance Explained in Y');
На практике, вероятно, было бы целесообразно проявлять большую осторожность при выборе количества компонентов. Перекрестная проверка, например, является широко используемым методом, который будет проиллюстрирован ниже в этом примере. На данный момент вышеприведенный график предполагает, что PLSR с двумя компонентами объясняет большую часть дисперсии наблюдаемого y. Вычислите подогнанные значения отклика для двухкомпонентной модели.
[Xloadings,Yloadings,Xscores,Yscores,betaPLS] = plsregress(X,y,2); yfitPLS = [ones(n,1) X]*betaPLS;
Затем поместите модель PCR с двумя основными компонентами. Первым шагом является выполнение анализа основных компонентов для X, с использованием pca и сохранение двух главных компонентов. Затем ПЦР представляет собой просто линейную регрессию переменной ответа на эти два компонента. Часто имеет смысл сначала нормализовать каждую переменную по ее стандартному отклонению, когда переменные имеют очень разные величины изменчивости, однако это не делается здесь.
[PCALoadings,PCAScores,PCAVar] = pca(X,'Economy',false);
betaPCR = regress(y-mean(y), PCAScores(:,1:2));
Чтобы облегчить интерпретацию результатов ПЦР в терминах исходных спектральных данных, преобразуйте их в коэффициенты регрессии для исходных незацентрированных переменных.
betaPCR = PCALoadings(:,1:2)*betaPCR; betaPCR = [mean(y) - mean(X)*betaPCR; betaPCR]; yfitPCR = [ones(n,1) X]*betaPCR;
График подобран по сравнению с наблюдаемым ответом для PLSR и PCR подходит.
plot(y,yfitPLS,'bo',y,yfitPCR,'r^'); xlabel('Observed Response'); ylabel('Fitted Response'); legend({'PLSR with 2 Components' 'PCR with 2 Components'}, ... 'location','NW');
В некотором смысле сравнение на графике выше не является справедливым - количество компонентов (два) было выбрано, глядя на то, насколько хорошо двухкомпонентная модель PLSR предсказала ответ, и нет причин, почему модель PCR должна быть ограничена тем же количеством компонентов. Однако при одинаковом количестве компонентов PLSR выполняет гораздо более эффективную работу по установке y. На самом деле, глядя на горизонтальный разброс подогнанных значений на графике выше, ПЦР с двумя компонентами едва ли лучше, чем использование постоянной модели. Значения r-квадрата из двух регрессий подтверждают это.
TSS = sum((y-mean(y)).^2); RSS_PLS = sum((y-yfitPLS).^2); rsquaredPLS = 1 - RSS_PLS/TSS
rsquaredPLS =
0.9466
RSS_PCR = sum((y-yfitPCR).^2); rsquaredPCR = 1 - RSS_PCR/TSS
rsquaredPCR =
0.1962
Другим способом сравнения прогнозирующей силы двух моделей является построение графика переменной ответа против двух предикторов в обоих случаях.
plot3(Xscores(:,1),Xscores(:,2),y-mean(y),'bo'); legend('PLSR'); grid on; view(-30,30);
Это немного трудно увидеть, не имея возможности интерактивно вращать фигуру, но график PLSR выше показывает точки, близко разбросанные вокруг плоскости. С другой стороны, график ПЦР ниже показывает облако точек с небольшим указанием линейной зависимости.
plot3(PCAScores(:,1),PCAScores(:,2),y-mean(y),'r^'); legend('PCR'); grid on; view(-30,30);
Обратите внимание, что в то время как два компонента PLS являются гораздо лучшими предикторами наблюдаемого y, следующий рисунок показывает, что они объясняют несколько меньшую дисперсию наблюдаемого X чем первые два основных компонента, используемых в PCR.
plot(1:10,100*cumsum(PLSPctVar(1,:)),'b-o',1:10, ... 100*cumsum(PCAVar(1:10))/sum(PCAVar(1:10)),'r-^'); xlabel('Number of Principal Components'); ylabel('Percent Variance Explained in X'); legend({'PLSR' 'PCR'},'location','SE');
Тот факт, что кривая PCR равномерно выше, говорит о том, почему PCR с двумя компонентами выполняет такую плохую работу по сравнению с PLSR при установке y. PCR конструирует компоненты, чтобы лучше объяснить X, и в результате эти первые два компонента игнорируют информацию в данных, которая важна для подгонки наблюдаемых y.
По мере добавления большего количества компонентов в PCR, он обязательно будет лучше подгонять исходные данные. y, просто потому, что в какой-то момент большая часть важной прогностической информации в X будет присутствовать в основных компонентах. Например, следующий рисунок показывает, что разница в остатках для двух методов гораздо менее драматична при использовании десяти компонентов, чем для двух компонентов.
yfitPLS10 = [ones(n,1) X]*betaPLS10; betaPCR10 = regress(y-mean(y), PCAScores(:,1:10)); betaPCR10 = PCALoadings(:,1:10)*betaPCR10; betaPCR10 = [mean(y) - mean(X)*betaPCR10; betaPCR10]; yfitPCR10 = [ones(n,1) X]*betaPCR10; plot(y,yfitPLS10,'bo',y,yfitPCR10,'r^'); xlabel('Observed Response'); ylabel('Fitted Response'); legend({'PLSR with 10 components' 'PCR with 10 Components'}, ... 'location','NW');
Обе модели подходят y довольно точно, хотя PLSR все же делает несколько более точную подгонку. Однако десять компонентов все еще являются произвольно выбранным числом для любой модели.
Часто полезно выбирать количество компонентов, чтобы минимизировать ожидаемую ошибку при прогнозировании ответа из будущих наблюдений за переменными предиктора. Простое использование большого количества компонентов поможет подгонять текущие наблюдаемые данные, но является стратегией, которая приводит к переоборудованию. Слишком хорошая подгонка текущих данных приводит к модели, которая плохо обобщается с другими данными и дает чрезмерно оптимистичную оценку ожидаемой ошибки.
Перекрестная проверка является более статистически обоснованным методом выбора количества компонентов в PLSR или PCR. Это позволяет избежать переоборудования данных, не прибегая к повторному использованию одних и тех же данных как для соответствия модели, так и для оценки ошибки прогнозирования. Таким образом, оценка ошибки предсказания не является оптимистично смещенной вниз.
plsregress имеет возможность оценить среднеквадратичную ошибку прогнозирования (MSEP) путем перекрестной проверки, в данном случае с использованием 10-кратного C-V.
[Xl,Yl,Xs,Ys,beta,pctVar,PLSmsep] = plsregress(X,y,10,'CV',10);
Для PCR, crossval в сочетании с простой функцией для вычисления суммы квадратичных ошибок для PCR может оценить MSEP, снова используя 10-кратную перекрестную проверку.
PCRmsep = sum(crossval(@pcrsse,X,y,'KFold',10),1) / n;
Кривая MSEP для PLSR показывает, что два или три компонента выполняют максимально хорошую работу. С другой стороны, ПЦР необходимо четыре компонента, чтобы получить одинаковую точность прогнозирования.
plot(0:10,PLSmsep(2,:),'b-o',0:10,PCRmsep,'r-^'); xlabel('Number of components'); ylabel('Estimated Mean Squared Prediction Error'); legend({'PLSR' 'PCR'},'location','NE');
Фактически, второй компонент в ПЦР увеличивает ошибку прогнозирования модели, предполагая, что комбинация переменных предиктора, содержащихся в этом компоненте, не сильно коррелирует с y. Опять же, это потому, что ПЦР конструирует компоненты, чтобы объяснить вариацию в X, не y.
Таким образом, если для PCR требуется четыре компонента, чтобы получить такую же точность прогнозирования, как PLSR с тремя компонентами, является ли модель PLSR более разборной? Это зависит от того, какой аспект модели вы рассматриваете.
Веса PLS представляют собой линейные комбинации исходных переменных, которые определяют компоненты PLS, т.е. они описывают, насколько сильно каждый компонент в PLSR зависит от исходных переменных и в каком направлении.
[Xl,Yl,Xs,Ys,beta,pctVar,mse,stats] = plsregress(X,y,3); plot(1:401,stats.W,'-'); xlabel('Variable'); ylabel('PLS Weight'); legend({'1st Component' '2nd Component' '3rd Component'}, ... 'location','NW');
Аналогично, загрузки PCA описывают, насколько сильно каждый компонент в PCR зависит от исходных переменных.
plot(1:401,PCALoadings(:,1:4),'-'); xlabel('Variable'); ylabel('PCA Loading'); legend({'1st Component' '2nd Component' '3rd Component' ... '4th Component'},'location','NW');
Для PLSR или PCR может быть, что каждому компоненту может быть дана физически значимая интерпретация путем проверки, какие переменные он весит наиболее сильно. Например, с помощью этих спектральных данных можно интерпретировать пики интенсивности в терминах соединений, присутствующих в бензине, и затем наблюдать, что веса для конкретного компонента выбирают небольшое количество этих соединений. С этой точки зрения, меньше компонентов проще интерпретировать, и поскольку PLSR часто требует меньше компонентов для адекватного прогнозирования ответа, это приводит к более экономным моделям.
С другой стороны, и PLSR, и PCR приводят к одному коэффициенту регрессии для каждой из исходных переменных предиктора плюс перехват. В этом смысле ни одна из них не является более разумной, поскольку независимо от того, сколько компонентов используется, обе модели зависят от всех предикторов. Более конкретно, для этих данных обеим моделям требуется 401 спектральное значение интенсивности, чтобы сделать прогноз.
Однако конечная цель может заключаться в уменьшении первоначального набора переменных до меньшего подмножества, все еще способного точно предсказывать отклик. Например, можно использовать веса PLS или нагрузки PCA для выбора только тех переменных, которые вносят наибольший вклад в каждый компонент. Как показано ранее, некоторые компоненты из подгонки модели ПЦР могут служить в первую очередь для описания изменения переменных предиктора и могут включать большие веса для переменных, которые не сильно коррелируют с ответом. Таким образом, ПЦР может привести к сохранению переменных, которые не нужны для прогнозирования.
Для данных, используемых в этом примере, разница в количестве компонентов, необходимых PLSR и PCR для точного прогнозирования, не велика, и веса PLS и нагрузки PCA, по-видимому, выбирают одни и те же переменные. Это может быть неверно для других данных.