[R-br] Como recriar um componente de PCA obtido com svyprcomp()

Leonardo Ferreira Fontenelle leonardof em leonardof.med.br
Domingo Abril 24 17:13:59 BRT 2016


Parece que não sou a única pessoa com dificuldade em gerar esse vetor
(quase) perfeitamente correlacionado com o primeiro vetor. O próprio
comando predict não consegue!
 
Retomando os exemplos anteriores:
 
pca1 <- prcomp(USArrests, scale = TRUE)
predpca1 <- predict(pca1, USArrests)
cor(pca1$x[, "PC1"], predpca1[, "PC1"])
# resultado: 1,00
 
library("survey")
data(api)
dclus2 <- svydesign(id = ~ dnum + snum, fpc = ~ fpc1 + fpc2, data
= apiclus2)
pca2 <- svyprcomp(~ api99 + api00 + ell, design = dclus2, scale = TRUE,
scores = TRUE)
predpca2 <- predict(pca2, dclus2$variables)
cor(pca2$x[, "PC1"], predpca2[, "PC1"])
# resultado: 0,506
 
Isso é um erro do  svyprcomp??
 
Leonardo Ferreira Fontenelle[1]
 
 
Em Dom 24 abr. 2016, às 16:11, Leonardo Ferreira Fontenelle escreveu:
> Em Dom 24 abr. 2016, às 00:16, Leonardo Ferreira Fontenelle escreveu:
>> Boa tarde a todos!
>>
>> Estou tentando escrever os resultados de uma análise de componentes
>> principais de forma que os leitores possam calcular o primeiro
>> componente a partir dos dados sem a necessidade de utilizar um
>> comando de análise de componentes principais propriamente dito. Minha
>> ideia é informar a codificação das variáveis, as cargas, os desvios-
>> padrão, e então orientar os leitores a multiplicar os valores das
>> variáveis pelas cargas e dividir pelos desvios-padrão e então somar
>> tudo. Minha expectativa é de que o resultado seja altamente
>> correlacionado àquele obtido por uma análise de componentes
>> principais propriamente dita. No entanto, como explico a seguir, não
>> estou conseguindo, e de alguma forma isso tem a ver com eu estar
>> usando svyprcomp() em vez de prcomp().
>>
>> Para começar, mostro uma situação em que isso funciona:
>>
>> pca1 <- prcomp(USArrests, scale = TRUE)
>> table1 <- data.frame(loadings = pca1$rotation[, "PC1"], scale =
>> pca1$scale, coef = NA_real_, row.names = row.names(pca1$rotation))
>> table1$coef <- table1$loadings / table1$scale
>> firstcomponent1 <- with(USArrests, Murder * table1["Murder", "coef"]
>> + Assault * table1["Assault", "coef"] + UrbanPop * table1["UrbanPop",
>> "coef"] + Rape * table1["Rape", "coef"])
>> cor(firstcomponent1, pca1$x[, "PC1"])
>> # resultado: 1,00
>>
>> Agora, uma situação onde isso não funciona:
>>
>> library("survey")
>> data(api)
>> dclus2 <- svydesign(id = ~ dnum + snum, fpc = ~ fpc1 + fpc2, data =
>> apiclus2)
>> pca2 <- svyprcomp(~ api99 + api00 + ell, design = dclus2, scale =
>> TRUE, scores = TRUE)
>> table2 <- data.frame(loadings = pca2$rotation[, "PC1"], scale =
>> pca2$scale, coef = NA_real_, row.names = row.names(pca2$rotation))
>> table2$coef <- table2$loadings / table2$scale
>> firstcomponent2 <- with(apiclus2, api99 * table2["api99", "coef"] +
>> api00 * table2["api00", "coef"] + ell * table2["ell", "coef"])
>> cor(firstcomponent2, pca2$x[, "PC1"])
>> # resultado: 0,506
>>
>> Como eu consigo cargas ("loadings") e desvios-padrão ("scale") que eu
>> possa usar de forma a conseguir calcular à mão algum vetor altamente
>> correlacionado com o primeiro componente do svyprcomp()?
>>
>> Grato!
>>
>> Leonardo Ferreira Fontenelle[2]
>>
>
> O problema tem a ver com a ponderação das observações. Não apenas essa
> parece ser a única diferença entre prcomp() e svyprcomp() (a função
> não leva em consideração correlações dentro dos conglomerados) mas
> também usar svyprcom() num delineamento gera componentes com uma
> correlação quase perfeita com os índices criados à mão como no e-mail
> anterior.
>
> Digitar "svyprcomp" na linha de comando do R, sem aspas ou parênteses,
> mostra o código-fonte. Ele me parece usar uma abordagem equivalente à
> da resposta http://stats.stackexchange.com/a/113488. Mesmo assim, não
> consigo entender de que forma posso criar um índice à mão que seja
> (quase) perfeitamente correlacionado com o primeiro componente.
>
> Leonardo Ferreira Fontenelle[3]
 

Links:

  1. http://lattes.cnpq.br/9234772336296638
  2. http://lattes.cnpq.br/9234772336296638
  3. http://lattes.cnpq.br/9234772336296638
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160424/47fec487/attachment.html>


Mais detalhes sobre a lista de discussão R-br