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

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

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

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

Заметьте, что в то время как два компонентов PLS являются намного лучшими предикторами наблюдаемой yследующие рисунки показывают, что они объясняют несколько меньшее отклонение наблюдаемых X чем первые два основных компонента, используемых в ПЦР.

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

Тот факт, что ПЦР-кривая равномерно выше, говорит о том, почему ПЦР с двумя компонентами делает такую плохую работу, относительно ПЛСР, в подборе кривой 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 или ПЦР. Он избегает избыточной подгонки данных, не используя те же данные повторно, чтобы соответствовать модели и оценить ошибку предсказания. Таким образом, оценка ошибки предсказания не является оптимистически смещенной вниз.

plsregress имеет опцию оценить среднюю квадратную ошибку предсказания (MSEP) путем перекрестной валидации, в этом случае используя 10-кратное C-V.

[Xl,Yl,Xs,Ys,beta,pctVar,PLSmsep] = plsregress(X,y,10,'CV',10);

Для ПЦР, crossval В сочетании с простой функцией для вычисления суммы квадратичных невязок для ПЦР, можно оценить 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');

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

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

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

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