Как выполнить ортогональную регрессию (наименьших квадратов) с помощью PCA?

Я всегда использую lm()в R для выполнения линейной регрессии $y$ на $x$ . Эта функция возвращает коэффициент $\beta$ такой, что

$$y = \beta x.$$

Сегодня я узнал об общих наименьших квадратах, и эту princomp()функцию (анализ основных компонентов, PCA) можно использовать для ее выполнения. Это должно быть хорошо для меня (точнее). Я сделал несколько тестов, используя princomp():

r <- princomp( ~ x + y)

Моя проблема: как интерпретировать ее результаты? Как я могу получить коэффициент регрессии? Под «коэффициентом» я подразумеваю число $\beta$ которое я должен использовать для умножения значения $x$ чтобы получить число, близкое к $y$ .

Обычные наименьшие квадраты против общих наименьших квадратов

Давайте сначала рассмотрим простейший случай только одной предикторной (независимой) переменной . Для простоты, пусть и x, и y центрированы, т. Е. Точка пересечения всегда равна нулю. Различие между стандартной регрессией OLS и «ортогональной» регрессией TLS ясно показано на этом (адаптированном мной) рисунке из самого популярного ответа в самой популярной теме на PCA:$x$$x$$y$

МНК припадки уравнение путем минимизации квадратов расстояний между наблюдаемыми значениями Y и прогнозируемых значений у . TLS соответствует тому же уравнению, минимизируя квадратные расстояния между точками ( x , y ) и их проекцией на линию. В этом простейшем случае линия TLS является просто первым основным компонентом двумерных данных. Чтобы найти β , проведите PCA по ( x , y ) точкам, т.е. построите$y=\beta x$$y$$\hat y$$(x,y)$$\beta$$(x,y)$ковариационную матрицу 2 × 2 Σ и найдите ее первый собственный вектор v$2\times 2$$\boldsymbol \Sigma$. ; тогда β = v y / v $\mathbf v = (v_x, v_y)$$\beta = v_y/v_x$

В Matlab:

 v = pca([x y]);    //# x and y are centered column vectors
 beta = v(2,1)/v(1,1);

В R:

 v <- prcomp(cbind(x,y))$rotation
 beta <- v[2,1]/v[1,1]

Кстати, это даст правильный наклон, даже если и y не были отцентрированы (потому что встроенные функции PCA автоматически выполняют центрирование). Чтобы восстановить перехват, вычислите β 0 = ˉ y - $x$$y$$\beta_0 = \bar y - \beta \bar x$ .

OLS против TLS, множественная регрессия

Для заданной зависимой переменной и множества независимых переменных x i (опять же, все по центру для простоты) регрессия соответствует уравнению y = β 1 x 1 + … + β p x p . МНК делает подгонку путем минимизации квадратов ошибок между наблюдаемыми значениями у и предсказанных значений у . TLS подгоняет, минимизируя квадрат расстояния между наблюдаемыми ( x , y ) ∈ R p + $y$$x_i$

$$y= \beta_1 x_1 + \ldots + \beta_p x_p.$$ $y$$\hat y$$(\mathbf x, y)\in\mathbb R^{p+1}$ точки и ближайшие точки на плоскости регрессии / гиперплоскости.

Обратите внимание, что больше нет «линии регрессии»! Вышеупомянутое уравнение определяет гиперплоскость : это двухмерная плоскость, если есть два предиктора, трехмерная гиперплоскость, если есть три предиктора, и т. Д. Таким образом, приведенное выше решение не работает: мы не можем получить решение TLS, взяв только первый ПК (который является линия). Тем не менее, решение может быть легко получено через PCA.

Как и раньше, PCA проводится по точкам. Это дает р + 1 собственных векторов в столбцах V . Первые р собственных векторов определяют р - мерную гиперплоскость H , что нам нужно; последний (число p + 1 ) собственный вектор v p + 1 ортогонален ему. Вопрос в том, как преобразовать основу H, заданную первым$(\mathbf x, y)$$p+1$$\mathbf V$$p$$p$$\mathcal H$$p+1$$\mathbf v_{p+1}$$\mathcal H$ собственных векторов в$p$$\boldsymbol \beta$ коэффициенты.

Заметим , что если положить для всех я ≠ K и только х К = 1 , то Y = р к , то есть вектор ( 0 , ... , 1 , ... , β к ) ∈ H лежит в гиперплоскости H , С другой стороны, мы знаем, что v p + 1 = ( v 1 , … , v p + $x_i=0$$i \ne k$$x_k=1$$\hat y=\beta_k$

$$(0,\ldots, 1, \ldots, \beta_k) \in \mathcal H$$ $\mathcal H$ ортогонально этому. Т.е. их скалярное произведение должно быть равно нулю: v k + β k v p + 1 = 0 ⇒ β k = - v k / v p + 1 $$\mathbf v_{p+1}=(v_1, \ldots, v_{p+1}) \:\bot\: \mathcal H$$ $$v_k + \beta_k v_{p+1}=0 \Rightarrow \beta_k = -v_k/v_{p+1}.$$

В Matlab:

 v = pca([X y]);    //# X is a centered n-times-p matrix, y is n-times-1 column vector
 beta = -v(1:end-1,end)/v(end,end);

В R:

 v <- prcomp(cbind(X,y))$rotation
 beta <- -v[-ncol(v),ncol(v)] / v[ncol(v),ncol(v)]

Опять же, это даст правильные наклоны, даже если и y не были отцентрированы (потому что встроенные функции PCA автоматически выполняют центрирование). Чтобы восстановить перехват, вычислите β 0 = ˉ y - ˉ x β .$x$$y$$\beta_0 = \bar y - \bar {\mathbf x} \boldsymbol \beta$

$x$$(x,y)$$v^{(1)}_y/v^{(1)}_x=-v^{(2)}_x/v^{(2)}_y$

Решение в закрытой форме для TLS

$\boldsymbol \beta$

$\mathbf X$$\mathbf y$$\mathbf v_{p+1}$$[\mathbf X\: \mathbf y]$$\sigma^2_{p+1}$$-\mathbf v_{p+1}/v_{p+1} = (\boldsymbol \beta\:\: -1)^\top$

$$\left(\begin{array}{c}\mathbf X^\top \mathbf X & \mathbf X^\top \mathbf y\\ \mathbf y^\top \mathbf X & \mathbf y^\top \mathbf y\end{array}\right) \left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right) = \sigma^2_{p+1}\left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right),$$ $$\boldsymbol \beta_\mathrm{TLS} = (\mathbf X^\top \mathbf X - \sigma^2_{p+1}\mathbf I)^{-1} \mathbf X^\top \mathbf y,$$ $$\boldsymbol \beta_\mathrm{OLS} = (\mathbf X^\top \mathbf X)^{-1} \mathbf X^\top \mathbf y.$$

Многофакторная множественная регрессия

Эту же формулу можно обобщить для многомерного случая, но даже для определения того, что делает многомерный TLS, потребуется некоторая алгебра. Смотрите Википедию на TLS . Многомерная регрессия OLS эквивалентна группе одномерных регрессий OLS для каждой зависимой переменной, но в случае TLS это не так.

Основываясь на найденной здесь наивной реализации GNU Octave , может сработать что-то вроде этого (крупица соли, уже поздно).

tls <- function(A, b){

  n <- ncol(A)
  C <- cbind(A, b)

  V <- svd(C)$v
  VAB <- V[1:n, (n+1):ncol(V)]
  VBB <- V[(n+1):nrow(V), (n+1):ncol(V)]
  return(-VAB/VBB)
}

princompвместо анализа регрессии общих наименьших квадратов выполняется анализ главных компонент . Насколько я знаю, нет ни функции R, ни пакета, который делает TLS; самое большее, это регрессия Деминга в MethComp .
Тем не менее, пожалуйста, воспринимайте это как предположение, что оно, скорее всего, не стоит того.