Загрузка данных

Загрузите набор данных, включающий спектральную интенсивность 60 выборок бензина в 401 длине волн и их оценок октана. Эти данные описаны в Kalivas, Джоне Х., "Два Набора данных Близких Инфракрасных Спектров", Хемометрики и Интеллектуальные Лабораторные Системы, 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 с десятью компонентами PLEASE и одним ответом.

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 функция и сохранение двух основных компонентов. PCR является затем только линейной регрессией переменной отклика на тех двух компонентах. Это часто целесообразно нормировать каждую переменную сначала на ее стандартное отклонение, когда переменные имеют совсем другие суммы изменчивости, однако, который не сделан здесь.

[PCALoadings,PCAScores,PCAVar] = pca(X,'Economy',false);
betaPCR = regress(y-mean(y), PCAScores(:,1:2));

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

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. На самом деле, смотря на горизонтальное рассеяние приспособленных значений в графике выше, PCR с двумя компонентами едва лучше, чем использование постоянной модели. r-squared значения от этих двух регрессий подтверждают это.

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 выше точек показов, тесно рассеянных о плоскости. С другой стороны, график PCR ниже показов облако точек с небольшой индикацией относительно линейного соотношения.

plot3(PCAScores(:,1),PCAScores(:,2),y-mean(y),'r^');
legend('PCR');
grid on; view(-30,30);

Заметьте, что, в то время как два компонента PLEASE являются намного лучшими предикторами наблюдаемого 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 указывает, что два или три компонента делают максимально хорошее задание. С другой стороны, PCR нужны четыре компонента, чтобы получить ту же точность предсказания.

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

На самом деле второй компонент в PCR увеличивает ошибку предсказания модели, предполагая, что комбинация переменных предикторов, содержавшихся в том компоненте, строго не коррелируется с y. Снова, поэтому PCR создает компоненты, чтобы объяснить изменение X, не y.

Бережливость модели

Таким образом, если PCR требует, чтобы четыре компонента получили ту же точность предсказания как PLSR с тремя компонентами, действительно ли модель PLSR более экономна? Это зависит от того, какой аспект модели вы рассматриваете.

Веса PLEASE являются линейными комбинациями исходных переменных, которые задают компоненты PLEASE, i.e., они описывают, как строго каждый компонент в 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, может случиться так, что каждому компоненту можно дать физически значимую интерпретацию путем осмотра, какие переменные это взвешивает наиболее в большой степени. Например, с этими спектральными данными может быть возможно интерпретировать peaks интенсивности в терминах составных объектов, существующих в бензине, и затем заметить, что веса для конкретного компонента выбирают небольшое количество тех составных объектов. С той точки зрения меньше компонентов более просто интерпретировать, и потому что PLSR часто требует, чтобы меньше компонентов предсказало ответ соответственно, это приводит к большему количеству экономных моделей.

С другой стороны, и PLSR и PCR приводят к одному коэффициенту регрессии для каждого из исходных переменных предикторов плюс точка пересечения. В этом смысле ни один не более экономен, потому что независимо от того, сколько компонентов используется, обе модели зависят от всех предикторов. Более конкретно, для этих данных, для обеих моделей нужно 401 спектральное значение интенсивности для того, чтобы сделать предсказание.

Однако май конечной цели, чтобы уменьшать исходный набор переменных к меньшему подмножеству, которое все еще в состоянии предсказать ответ точно. Например, может быть возможно использовать веса PLEASE или загрузки PCA, чтобы выбрать только те переменные, которые способствуют больше всего каждому компоненту. Как показано ранее некоторые компоненты от подгонки модели PCR могут служить, в основном, чтобы описать изменение переменных предикторов и могут включать большие веса для переменных, которые строго не коррелируются с ответом. Таким образом PCR может привести к сдерживающим переменным, которые являются ненужными для предсказания.

Для данных, используемых в этом примере, различие в количестве компонентов, необходимых PLSR и PCR для точного предсказания, не является большим, и веса PLEASE, и загрузки PCA, кажется, выбирают те же переменные. Это не может быть верно для других данных.