Regressão parcial de mínimos quadrados em R: por que o PLS em dados padronizados não é equivalente a maximizar a correlação?

12

Eu sou muito novo em mínimos quadrados parciais (PLS) e tento entender a saída da função R plsr()no plspacote. Vamos simular dados e executar o PLS:

library(pls)
n <- 50
x1 <- rnorm(n); xx1 <- scale(x1) 
x2 <- rnorm(n); xx2 <- scale(x2)
y <- x1 + x2 + rnorm(n,0,0.1); yy <- scale(y)
p <- plsr(yy ~ xx1+xx2, ncomp=1)

Eu estava esperando que os seguintes números a e b

> ( w <- loading.weights(p) )

Loadings:
    Comp 1
xx1 0.723 
xx2 0.690 

               Comp 1
SS loadings       1.0
Proportion Var    0.5
> a <- w["xx1",]
> b <- w["xx2",]
> a^2+b^2
[1] 1

são calculados para maximizar

> cor(y, a*xx1+b*xx2)
          [,1]
[1,] 0.9981291

mas este não é exatamente o caso:

> f <- function(ab){
+ a <- ab[1]; b <- ab[2]
+ cor(y, a*xx1+b*xx2)
+ }
> optim(c(0.7,0.6), f, control=list(fnscale=-1))
$par
[1] 0.7128259 0.6672870

$value
[1] 0.9981618

É um erro numérico, ou eu não compreendem a natureza de a e b ?

Eu também gostaria de saber quais são esses coeficientes:

> p$coef
, , 1 comps

           yy
xx1 0.6672848
xx2 0.6368604 

EDIT : Agora eu vejo o que p$coefé:

> x <- a*xx1+b*xx2
> coef(lm(yy~0+x))
        x 
0.9224208 
> coef(lm(yy~0+x))*a
        x 
0.6672848 
> coef(lm(yy~0+x))*b
        x 
0.6368604 

Então eu acho que estou certo sobre a natureza de ea .b

EDIT: Tendo em vista os comentários feitos por @chl, sinto que minha pergunta não é clara o suficiente, então deixe-me fornecer mais detalhes. Na minha exemplo, existe um vector de de respostas e uma matriz de duas colunas X de preditores e uso a versão normalizada ~ Y de Y e a versão normalizada ~ X de X (centrado e dividido pelo desvio padrão). A definição do primeiro componente PLS t 1 é t 1 = a ˜ X 1 + b ˜ X 2 com a e bYXY~YX~Xt1t1=aX~1+bX~2abescolhido de modo a ter um valor máximo do produto interno . Portanto, é equivalente a maximizar a correlação entre t 1 e Yt1,Y~t1Y , não é?

Stéphane Laurent
fonte
2
A regressão PLS maximiza a pontuação dos fatores (que são computados como produto de dados brutos com vetor (es) de carga) covariância , não correlação (como é feito na Análise de correlação canônica). Há uma boa visão geral do plspacote e da regressão PLS neste documento JSS .
chl
1
Como todos os vetores são centralizados e normalizados, covariância é correlação, não é? Desculpe, mas o artigo JSS é muito técnico para iniciantes.
Stéphane Laurent
Em geral, há um processo de deflação assimétrico (resultante da regressão da combinação linear de um bloco na combinação linear do outro) que complica um pouco as coisas. Forneci uma imagem esquemática nesta resposta . Hervé Abdi deu uma visão geral da regressão do PLS, e o Método de Levantamento dos Mínimos Quadrados Parciais (PLS) de Wegelin também é bastante útil. Neste ponto, eu provavelmente deveria converter todos esses comentários para uma resposta ...
chl
Na minha exemplo, existe um vector de de respostas e uma matriz de duas colunas X de preditores e uso a versão normalizada ~ Y de Y e a versão normalizada ~ X de X (centrado e dividido pelo desvio padrão). Meu definição dos primeiros PLS componente t 1 é t 1 = um ~ X 1 + b ~ X 2 com um e b escolhida para ter um valor máximo do produto escalar t 1 ,YXY~YX~Xt1t1=aX~1+bX~2ab. Não é a boa definição? t1,Y~
Stéphane Laurent
Desculpe, @ Stéphane, porque meus comentários acima não levaram em conta o fato de que você solicitou apenas um componente (para que a deflação não desempenhe um papel crítico aqui). No entanto, parece que a sua função de otimização não impõe vetores de pesos norma unidade, de modo que no final. (aliás, lhe dará mais informações sobre o que são esses 'coeficientes', mas você já descobriu que se aparentemente.)a2+b21?coef.mvr
chl

Respostas:

17

A regressão PLS depende de algoritmos iterativos (por exemplo, NIPALS, SIMPLS). Sua descrição das idéias principais está correta: buscamos um (PLS1, uma variável de resposta / múltiplos preditores) ou dois (PLS2, com modos diferentes, variáveis ​​de resposta múltipla / múltiplos preditores) vetor (es) de pesos, (e v ) digamos, para formar combinações lineares da (s) variável (s) original (ais), de modo que a covariância entre Xu e Y (Yv, para PLS2) seja máxima. Vamos nos concentrar na extração do primeiro par de pesos associado ao primeiro componente. Formalmente, o critério para otimizar lê max cov ( X u , Yuv No seu caso, Y é univariada, portanto, eleva-se a maximizar cov ( X u , y ) Var ( X u ) 1 / 2 x CR ( X u , y ) × Var ( Y ) 1 / 2

maxcov(Xvocê,Yv).(1)
Y Desde Var ( y ) não depende L , temos que maximizar Var ( X u ) 1 / 2 x CR ( X u , y ) . Vamos considerar, onde os dados são padronizados individualmente (inicialmente cometi o erro de dimensionar sua combinação linear em vez de x 1 e x 2 separadamente!), Para que Var ( x Var (
cov(Xvocê,y)Var(Xvocê)1/2×cor(Xvocê,y)×Var(y)1/2,st.__você__=1
Var(y)vocêVar(Xvocê)1/2×cor(Xvocê,y)X=[x_1;x_2]x1x2 ; no entanto, Var ( X u ) 1 e depende de uVar(x1)=Var(x2)=1Var(Xvocê)1você . Em conclusão, maximizar a correlação entre o componente latente e a variável de resposta não produzirá os mesmos resultados .

Eu deveria agradecer Arthur Tenenhaus que me indicou a direção certa.

O uso de vetores de peso unitário não é restritivo e alguns pacotes ( pls. regressionem plsgenomics , com base no código do pacote anterior de Wehrens pls.pcr) retornarão vetores de peso não padronizados (mas com componentes latentes ainda da norma 1), se solicitado. Mas a maioria dos pacotes PLS retornará u padronizadosvocê , incluindo o que você usou, principalmente aqueles que implementam o algoritmo SIMPLS ou NIPALS; Encontrei uma boa visão geral de ambas as abordagens na apresentação de Barry M. Wise, Propriedades da regressão de mínimos quadrados parciais (PLS) e diferenças entre algoritmos , mas a quimiometria.a vinheta também oferece uma boa discussão (págs. 26-29). Também é de particular importância o fato de que a maioria das rotinas PLS (pelo menos a que conheço em R) presume que você forneça variáveis ​​não padronizadas porque a centralização e / ou o dimensionamento são tratados internamente (isso é particularmente importante ao realizar a validação cruzada, por exemplo )

Dada a restrição , o vetor u é u = X yuu=1u

u=XyXy.

Usando um pouco de simulação, ele pode ser obtido da seguinte maneira:

set.seed(101)
X <- replicate(2, rnorm(100))
y <- 0.6*X[,1] + 0.7*X[,2] + rnorm(100)
X <- apply(X, 2, scale)
y <- scale(y)

# NIPALS (PLS1)
u <- crossprod(X, y)
u <- u/drop(sqrt(crossprod(u)))         # X weights
t  <- X%*%u
p <- crossprod(X, t)/drop(crossprod(t)) # X loadings

Você pode comparar os resultados acima ( u=[0.5792043;0.8151824]em particular) com o que os pacotes R dariam. Por exemplo, usando NIPALS do pacote quimiométrico (outra implementação que eu sei que está disponível no pacote mixOmics ), obteríamos :

library(chemometrics)
pls1_nipals(X, y, 1)$W  # X weights [0.5792043;0.8151824]
pls1_nipals(X, y, 1)$P  # X loadings

Resultados semelhantes seriam obtidos com plsro algoritmo PLS do kernel padrão:

> library(pls)
> as.numeric(loading.weights(plsr(y ~ X, ncomp=1)))
[1] 0.5792043 0.8151824

u tem comprimento 1.

Desde que você altere sua função para otimizar para uma que leia

f <- function(u) cov(y, X%*%(u/sqrt(crossprod(u))))

e normalizar udepois (u <- u/sqrt(crossprod(u)) ), você deve estar mais perto da solução acima.

maxuXYv,
uXY
svd(crossprod(X, y))$u

No caso mais geral (PLS2), uma maneira de resumir o exposto acima é dizer que os primeiros vetores canônicos do PLS são a melhor aproximação da matriz de covariância de X e Y em ambas as direções.

Referências

  1. Tenenhaus, M (1999). L'approche PLS . Revue de Statistique Appliquée , 47 (2), 5-40.
  2. ter Braak, CJF e de Jong, S (1993). A função objetiva da regressão parcial de mínimos quadrados . Journal of Chemometrics , 12, 41-54.
  3. Abdi, H (2010). Regressão e projeção de mínimos quadrados parciais na regressão de estrutura latente (regressão PLS) . Wiley Interdisciplinary Reviews: Computational Statistics , 2, 97-106.
  4. Boulesteix, AL e Strimmer, K (2007). Mínimos quadrados parciais: uma ferramenta versátil para a análise de dados genômicos de alta dimensão . Briefings in Bioinformtics , 8 (1), 32-44.
chl
fonte
Obrigado chl. Vou ler a sua resposta sempre que possível (e, certamente, upvote e clique na marca de seleção!)
Stéphane Laurent
Acabei de ler sua resposta - parabéns e muito obrigado.
Stéphane Laurent