Como executar a regressão ortogonal (total de mínimos quadrados) via PCA?

29

Eu sempre uso lm()em R para executar regressão linear de y em . Essa função retorna um coeficiente tal quexβ

y=βx.

Hoje eu aprendi sobre o total de mínimos quadrados e essa princomp()função (análise de componentes principais, PCA) pode ser usada para realizá-lo. Deve ser bom para mim (mais preciso). Eu fiz alguns testes usando princomp(), como:

r <- princomp( ~ x + y)

Meu problema é: como interpretar seus resultados? Como posso obter o coeficiente de regressão? Por "coeficiente", quero dizer o número que eu tenho que usar para multiplicar o valor para fornecer um número próximo de .x yβxy

Dail
fonte
Um momento pessoal, estou um pouco confuso. veja: zoonek2.free.fr/UNIX/48_R/09.html Isso é chamado de PCA (Análise de componentes principais, também conhecida como "regressão ortogonal" ou "somas perpendiculares de quadrados" ou "total de mínimos quadrados"), então acho que estamos falando sobre TLS com princomp () Não?
Dail
Não; essas são duas coisas diferentes, consulte o artigo da wikipedia sobre o PCA. O fato de ser usado aqui é um hack (eu não sei o quão exato, mas vou verificar); é por isso que a extração complexa de coeficientes.
1
Uma pergunta relacionada: stats.stackexchange.com/questions/2691/… e uma postagem no blog são referenciadas por uma das respostas: cerebralmastication.com/2010/09/…
Jonathan

Respostas:

48

Mínimos quadrados comuns vs. mínimos quadrados totais

Vamos primeiro considerar o caso mais simples de apenas uma variável preditora (independente) . Para simplificar, deixe x e y centralizados, ou seja, a interceptação é sempre zero. A diferença entre a regressão OLS padrão e a regressão TLS "ortogonal" é mostrada claramente nesta figura (adaptada por mim) da resposta mais popular no segmento mais popular no PCA:xxy

OLS vs TLS

OLS ajusta a equação y=βx minimizando distâncias quadradas entre os valores observados y e valores preditos y . TLS se encaixa na mesma equação, minimizando distâncias ao quadrado entre ( x , yy^ e sua projeção na linha. Neste caso mais simples, a linha TLS é simplesmente o primeiro componente principal dos dados 2D. Para encontrar β , do APC em ( x , y ) pontos, isto é, a construção de 2 × 2 covariância matriz Σ e encontrar o seu primeiro vector próprio v =(x,y)β(x,y)2×2Σ ; então β = v y / v x .v=(vx,vy)β=vy/vx

No Matlab:

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

Em R:

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

A propósito, esta vai produzir inclinação correcta, mesmo que x e y não foram centrado (porque funções internas PCA executar automaticamente centragem). Para recuperar a interceptação, calcule .β0=y¯βx¯

OLS vs. TLS, regressão múltipla

Dada uma variável dependente e muitas variáveis ​​independentes x i (novamente, todas centradas na simplicidade), a regressão se ajusta a uma equação y = β 1 x 1 + + β p x p . O OLS faz o ajuste minimizando os erros ao quadrado entre os valores observados de y e os valores previstosyxi

y=β1x1++βpxp.
y . O TLS faz o ajuste minimizando as distâncias ao quadrado entre as observações(x,y)Rp+1y^(x,y)Rp+1 pontos e os pontos mais próximos no plano de regressão / hiperplano.

Observe que não há mais "linha de regressão"! A equação acima especifica um hiperplano : é um plano 2D se houver dois preditores, hiperplano 3D se houver três preditores, etc. Portanto, a solução acima não funciona: não podemos obter a solução TLS usando apenas o primeiro PC (que é uma linha). Ainda, a solução pode ser facilmente obtida via PCA.

Como antes, o PCA é executado em pontos . Isto produz p + 1 vectores próprios em colunas de V . As primeiras p vectores próprios definir uma p -dimensional hiperplana H que é necessário; o último (número p + 1 ) autovetor(x,y)p+1VppHp+1 é ortogonal a ele. A questão é como transformar a base de H dado pelos primeirospvectores próprios para ospcoeficientes.vp+1Hpβ

Observe que, se definir para todos os i k e só x k = 1 , então y = p k , ou seja, o vector ( 0 , ... , 1 , ... , β k ) H reside no hiperplana H . Por outro lado, sabemos que vxi=0ikxk=1y^=βk

(0,,1,,βk)H
H é ortogonal a ele. Ou seja, seu produto escalar deve ser zero: v k + β k v p + 1 = 0 β k = - v k / v p + 1 .
vp+1=(v1,,vp+1)H
vk+βkvp+1=0βk=vk/vp+1.

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

Em R:

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

Mais uma vez, isto irá proporcionar pistas correctas, mesmo que e y não foram centrado (porque funções internas PCA executar automaticamente centragem). Para recuperar a interceptação, calcule β 0 = ˉ y - ˉ x β .xyβ0=y¯x¯β

Como verificação de integridade, observe que esta solução coincide com a anterior no caso de apenas um único preditor . De fato, o espaço ( x , y ) é 2D e, portanto, dado que o primeiro vetor próprio PCA é ortogonal ao segundo (último), v ( 1 ) y / v ( 1 ) x = - v ( 2 ) x / v ( 2 ) y .x(x,y)vy(1)/vx(1)=vx(2)/vy(2)

Solução de formulário fechado para TLS

Surpreendentemente, verifica-se que existe uma equação de forma fechada para . O argumento abaixo é retirado do livro de Sabine van Huffel "O total de mínimos quadrados" (seção 2.3.2).β

Seja e y as matrizes de dados centralizadas. O último vetor próprio PCAXy é um vetor próprio da matriz de covariância de[ Xvp+1 com um valor próprio σ 2 p + 1 . Se é um vetor próprio, então o é - v p +[Xy]σp+12 . Escrever a equação eigenvector: ( XX Xy yX yy ) ( β - 1 ) = σ 2 p + 1 ( β - 1 ) , e calcular o produto à esquerda, nós imediatamente obter esse β T L S = ( XX - σvp+1/vp+1=(β1)

(XXXyyXyy)(β1)=σp+12(β1),
que lembra fortemente a expressão familiar de OLS β O L S =( XX ) - 1 Xy .
βTLS=(XXσp+12I)1Xy,
βOLS=(XX)1Xy.

Regressão múltipla multivariada

A mesma fórmula pode ser generalizada para o caso multivariado, mas mesmo para definir o que o TLS multivariado faz, exigiria alguma álgebra. Veja a Wikipedia sobre TLS . A regressão OLS multivariada é equivalente a várias regressões OLS univariadas para cada variável dependente, mas, no caso do TLS, não é assim.

ameba diz Restabelecer Monica
fonte
1
Eu não conheço R, mas ainda queria fornecer trechos de R para referência futura. Há muitas pessoas aqui com proficiência em R. Por favor, sinta-se à vontade para editar meus trechos, se necessário! Obrigado.
Ameba diz Reinstate Monica
Bom post, mas se posso perguntar o que garante o fato de que o vetor está no hiperplano? (0,,1,,βk)
JohnK
xixk=1y=βjxjy=βk1=βk(0,,1,βk)y=βjxj
Ameba diz Reinstate Monica
Eu pareço ter interpretado mal essa parte, mas agora está claro. Obrigado pelo esclarecimento também.
JohnK
2
Em R, você pode preferir "eigen (cov (cbind (x, y))) $ vectors" em vez de "prcomp (cbind (x, y)) $ rotação" porque o primeiro é muito mais rápido para vetores maiores.
Thomas Browne
9

Com base na ingênua implementação GNU Octave encontrada aqui , algo como isto pode (grão de sal, é tarde) funcionar.

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)
}
cashoes
fonte
4

princompestá executando a análise de componentes principais em vez da regressão total de mínimos quadrados. Até onde eu sei, não há função R nem pacote que faça TLS; no máximo, há regressão de Deming no MethComp .
No entanto, trate isso como uma sugestão de que provavelmente não vale a pena.


fonte
Eu pensei que Deming no pacote MethComp fosse TLS - qual é a diferença?
mark999
Você deve fornecer a taxa de erros em xey; TLS puro otimiza isso.