Документация

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

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