Como recriar um componente de PCA obtido com svyprcomp()
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[1] Links: 1. http://lattes.cnpq.br/9234772336296638
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[1]
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[2] Links: 1. http://lattes.cnpq.br/9234772336296638 2. http://lattes.cnpq.br/9234772336296638
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
Leonardo, Você já está no terceiro post sobre o seu problema sem que ninguém tenha feito uma intervenção, é o que a gente chama de « depuração confessional » (veja https://books.google.com.br/books?id=4miO-X83hmUC&pg=PA255&lpg=PA255&dq=C+pr... )!! Sendo específico com relação ao que você reporta, acho que você deve tentar entender se o problema é com os dados ou com os métodos: Escolha um dos conjuntos de dados, ou o "USArrest" ou o "api" e aplique as duas formas de fazer a ACP, analise os resultados. Por último, e imagino que não seja importante para você desvendar o mistério que lhe apareceu, coloco a seguinte dúvida, qual a vantagem de fazer toda essa ginástica se para gerar o tal resultado que não necessitaria do comando para " análise de componentes principais propriamente dito" você usa coeficientes, loadings, etc. obtidos desse comando? HTH 2016-04-24 17:13 GMT-03:00 Leonardo Ferreira Fontenelle < leonardof@leonardof.med.br>:
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 <http://lattes.cnpq.br/9234772336296638>
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 <http://lattes.cnpq.br/9234772336296638>
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 <http://lattes.cnpq.br/9234772336296638>
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
O problema não é específico dos dados. Usei o "api" para facilitar a reprodução, mas esbarrei no problema ao analisar o banco de dados de um inquérito domiciliar de abrangência nacional. Minha intenção era criar um indicador econômico baseado em bens, que as pessoas poderiam utilizar em seus próprios estudos. Dizer "pega essas variáveis e joga no PCA" poderia até ser uma alternativa tecnicamente correta, mas restringiria o uso do indicador a quem sabe o que é isso. Dessa forma, seguindo o passo de outras iniciativas semelhantes, eu gostaria de dizer "pega as variáveis, multiplica por esse número e você tem o indicador XYZ". Outro motivo para eu poder dizer que o centro da questão é a ferramente é que eu rodei a PCA (com os pesos amostrais) no Stata, e os resultados foram semelhantes aos que eu consigo no R com o predict() ou com o "with(apiclus2, api99 * table2["api99", "coef"] + ... ", mas diferentes dos que eu consigo no R com o pca2$x[, "PC1"]. Dessa forma, embora a resposta possa admitidamente ser mas criativa do que isso, estou entre as possibilidades de que (a) os componentes principais do svyprcomp têm alguma interpretação especial que me escapa até agora, ou (b) o svyprcomp tem algum erro, que deveria ser identificado e reparado. Infelizmente, na estatística eu sou usário, então não consigo abrir o código do svyprcomp e dizer se a forma como ele lida com os pesos amostrais (ou a SVD) está correta ou não. Grato pela atenção, Leonardo Ferreira Fontenelle[1] Em Ter 26 abr. 2016, às 18:18, Cesar Rabak escreveu:
Leonardo,
Você já está no terceiro post sobre o seu problema sem que ninguém tenha feito uma intervenção, é o que a gente chama de « depuração confessional » (veja https://books.google.com.br/books?id=4miO-X83hmUC&pg=PA255&lpg=PA255&dq=C+pr...
Sendo específico com relação ao que você reporta, acho que você deve tentar entender se o problema é com os dados ou com os métodos:
Escolha um dos conjuntos de dados, ou o "USArrest" ou o "api" e aplique as duas formas de fazer a ACP, analise os resultados.
Por último, e imagino que não seja importante para você desvendar o mistério que lhe apareceu, coloco a seguinte dúvida, qual a vantagem de fazer toda essa ginástica se para gerar o tal resultado que não necessitaria do comando para " análise de componentes principais propriamente dito" você usa coeficientes, loadings, etc. obtidos desse comando?
HTH
2016-04-24 17:13 GMT-03:00 Leonardo Ferreira Fontenelle <leonardof@leonardof.med.br>:
__ 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[2]
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[3]
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[4]
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Links: 1. http://lattes.cnpq.br/9234772336296638 2. http://lattes.cnpq.br/9234772336296638 3. http://lattes.cnpq.br/9234772336296638 4. http://lattes.cnpq.br/9234772336296638
Caro Leonardo, Seu post enseja vários assuntos: A respeito de « criar um indicador econômico baseado em bens », você já considerou usar o CCEB da ABEP http://www.abep.org/criterio-brasil? "pega essas variáveis e joga no PCA" vai ser sempre diferente de "pega as variáveis, multiplica por esse número e você tem o indicador XYZ" porque no primeiro caso um procedimento matemático criado para encontrar as dimensões com maior variabilidade nos dados seria aplicado enquanto no segundo você estaria criando um indicador baseado em uma média ponderada (no mínimo) ou uma combinação linear das variáveis com coeficientes fixos a priori. A diferença entre os dois procedimentos está, como você intui, no uso diferente dos pesos e no fato de o projeto do inquérito para os dados do exemplo usarem aglomerados (/clusters/), no caso deste exemplo 15 deles, e tamanho das populações de 757 casos (essa informação é usada para corrigir a variância calculada [a famosa correção para população finita]). Por isso, minha humilde insistência que você não compare duas coisas diferentes, pois o exemplo do USArrests com procomp não tem nenhuma dessas sofisticações. Não entendo esta parte: «. . . os resultados foram semelhantes aos que eu consigo no R com o predict() ou com o "with(apiclus2, api99 * table2["api99", "coef"] + ... ", mas diferentes dos que eu consigo no R com o pca2$x[, "PC1"]. » Como você calculou os "coef"'s em cada SW? Para dizer se a svyprcomp pode ter algum erro ou não, você precisa criar um design que seja apropriado para poder usar os dados os USArrests na svyprcomp e conseguir chegar nos resultados da prcomp do pacote "base" do R, que indicaria que não há erro, ou se o condicionamento dos dados estiver correto para svyprcomp, então discutir a discrepância com os autores do pacote survey. [] 2016-04-26 21:56 GMT-03:00 Leonardo Ferreira Fontenelle < leonardof@leonardof.med.br>:
O problema não é específico dos dados. Usei o "api" para facilitar a reprodução, mas esbarrei no problema ao analisar o banco de dados de um inquérito domiciliar de abrangência nacional. Minha intenção era criar um indicador econômico baseado em bens, que as pessoas poderiam utilizar em seus próprios estudos. Dizer "pega essas variáveis e joga no PCA" poderia até ser uma alternativa tecnicamente correta, mas restringiria o uso do indicador a quem sabe o que é isso. Dessa forma, seguindo o passo de outras iniciativas semelhantes, eu gostaria de dizer "pega as variáveis, multiplica por esse número e você tem o indicador XYZ".
Outro motivo para eu poder dizer que o centro da questão é a ferramente é que eu rodei a PCA (com os pesos amostrais) no Stata, e os resultados foram semelhantes aos que eu consigo no R com o predict() ou com o "with(apiclus2, api99 * table2["api99", "coef"] + ... ", mas diferentes dos que eu consigo no R com o pca2$x[, "PC1"].
Dessa forma, embora a resposta possa admitidamente ser mas criativa do que isso, estou entre as possibilidades de que (a) os componentes principais do svyprcomp têm alguma interpretação especial que me escapa até agora, ou (b) o svyprcomp tem algum erro, que deveria ser identificado e reparado. Infelizmente, na estatística eu sou usário, então não consigo abrir o código do svyprcomp e dizer se a forma como ele lida com os pesos amostrais (ou a SVD) está correta ou não.
Grato pela atenção,
Leonardo Ferreira Fontenelle <http://lattes.cnpq.br/9234772336296638>
Em Ter 26 abr. 2016, às 18:18, Cesar Rabak escreveu:
Leonardo,
Você já está no terceiro post sobre o seu problema sem que ninguém tenha feito uma intervenção, é o que a gente chama de « depuração confessional » (veja https://books.google.com.br/books?id=4miO-X83hmUC&pg=PA255&lpg=PA255&dq=C+pr... )!!
Sendo específico com relação ao que você reporta, acho que você deve tentar entender se o problema é com os dados ou com os métodos:
Escolha um dos conjuntos de dados, ou o "USArrest" ou o "api" e aplique as duas formas de fazer a ACP, analise os resultados.
Por último, e imagino que não seja importante para você desvendar o mistério que lhe apareceu, coloco a seguinte dúvida, qual a vantagem de fazer toda essa ginástica se para gerar o tal resultado que não necessitaria do comando para " análise de componentes principais propriamente dito" você usa coeficientes, loadings, etc. obtidos desse comando?
HTH
2016-04-24 17:13 GMT-03:00 Leonardo Ferreira Fontenelle < leonardof@leonardof.med.br>:
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 <http://lattes.cnpq.br/9234772336296638>
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 <http://lattes.cnpq.br/9234772336296638>
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 <http://lattes.cnpq.br/9234772336296638>
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
*_______________________________________________* R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Em Ter 26 abr. 2016, às 22:52, Cesar Rabak escreveu:
Caro Leonardo,
Seu post enseja vários assuntos:
A respeito de « criar um indicador econômico baseado em bens », você já considerou usar o CCEB da ABEP http://www.abep.org/criterio-brasil?
Obrigado por estar ajudando, Cesar. Com relação ao CCEB, já até usei em estudo anterior (http://dx.doi.org/10.5712/rbmfc7(25)482[1]), mas neste caso ele não me atende. Além de eles serem ajustados a partir de levantamentos de dados restritos aos grandes centros urbanos, todas as versões usam um ou mais dados que estão ausentes no inquérito em questão. Conheço outro indicador econômico baseado em bens, mas ele foi criado para domicílios urbanos e também tem itens que não estão disponíveis no banco de dados.
"pega essas variáveis e joga no PCA" vai ser sempre diferente de "pega as variáveis, multiplica por esse número e você tem o indicador XYZ" porque no primeiro caso um procedimento matemático criado para encontrar as dimensões com maior variabilidade nos dados seria aplicado enquanto no segundo você estaria criando um indicador baseado em uma média ponderada (no mínimo) ou uma combinação linear das variáveis com coeficientes fixos a priori.
O "multiplica esse número" seria derivado da análise de componentes principais. Se fosse para criar um indicador sem PCA, eu não estaria quebrando a cabeça com essa questão :)
A diferença entre os dois procedimentos está, como você intui, no uso diferente dos pesos e no fato de o projeto do inquérito para os dados do exemplo usarem aglomerados (/clusters/), no caso deste exemplo 15 deles, e tamanho das populações de 757 casos (essa informação é usada para corrigir a variância calculada [a famosa correção para população finita]).
Por isso, minha humilde insistência que você não compare duas coisas diferentes, pois o exemplo do USArrests com procomp não tem nenhuma dessas sofisticações.
Olhando o código do svyprcomp(), ele não usa nem correção para população finita, nem conglomerados, apenas a ponderação das observações. Naturalmente, fazer PCA com e sem ponderação das observações tem que dar resultados diferentes, ainda mais em bancos de dados diferentes. O meu uso para o USArrests e o prcomp() foi mostrar que, nesse cenário mais usual, não existe essa discrepância entre métodos diferentes de obter o componente.
Não entendo esta parte: «. . . os resultados foram semelhantes aos que eu consigo no R com o predict() ou com o "with(apiclus2, api99 * table2["api99", "coef"] + ... ", mas diferentes dos que eu consigo no R com o pca2$x[, "PC1"]. »
Como você calculou os "coef"'s em cada SW?
Os coeficientes em questão são a divisão das cargas pelas escalas. Isso está ilustrado no código em meu primeiro e-mail, mas repito abaixo: 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 Existe uma forma mais elegante de obter os coeficientes, sem data.frame, mas só pensei nela depois de postar o primeiro e-mail. Com relação à correlação, eu até tenho no meu .Rprofile um "svycor" que calculava correlação a partir do svyvar[2] (), mas não achei pertinente usá-lo no código postado. Desde então, seguindo a "depuração confessional", encontrei uma alternativa mais facilmente reprodutível:
cov.wt(cbind(caseiro = firstcomponent2, original = pca2$x[, "PC1"]), wt = weights(dclus2), cor = TRUE)$cor # caseiro original # caseiro 1.0000000 0.4049466 # original 0.4049466 1.0000000
Naturalmente, esse comando não leva em consideração os conglomerados e FPC, mas ( como eu já disse) o próprio svyprcomp() também não.
Para dizer se a svyprcomp pode ter algum erro ou não, você precisa criar um design que seja apropriado para poder usar os dados os USArrests na svyprcomp e conseguir chegar nos resultados da prcomp do pacote "base" do R, que indicaria que não há erro, ou se o condicionamento dos dados estiver correto para svyprcomp, então discutir a discrepância com os autores do pacote survey.
Criando um objeto survey.design sem pesos amostrais nem conglomerados, o svyprcomp() se comporta igual ao prcomp():
d <- svydesign(ids = ~ 1, data = USArrests) Warning message: In svydesign.default(ids = ~1, data = USArrests) : No weights or probabilities supplied, assuming equal probability svypc <- svyprcomp(~ Murder + Assault + UrbanPop + Rape, d, scale. = TRUE, scores = TRUE) pca1 <- prcomp(USArrests, scale = TRUE) svypc Standard deviations: [1] 1.5748783 0.9948694 0.5971291 0.4164494
Rotation: PC1 PC2 PC3 PC4 Murder -0.5358995 0.4181809 -0.3412327 0.64922780 Assault -0.5831836 0.1879856 -0.2681484 -0.74340748 UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773 Rape -0.5434321 -0.1673186 0.8177779 0.08902432
pca1 Standard deviations: [1] 1.5748783 0.9948694 0.5971291 0.4164494
Rotation: PC1 PC2 PC3 PC4 Murder -0.5358995 0.4181809 -0.3412327 0.64922780 Assault -0.5831836 0.1879856 -0.2681484 -0.74340748 UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773 Rape -0.5434321 -0.1673186 0.8177779 0.08902432
cor(svypc$x[, "PC1"], predict(svypc, USArrests)[, "PC1"]) [1] 1
Não sei se ajuda na linha de raciocínio, mas também inverti a situação e apliquei o prcomp() ao api: dclus2 <- svydesign(id = ~ dnum + snum, fpc = ~ fpc1 + fpc2, data = apiclus2) svypc <- svyprcomp(~ api99 + api00 + ell, design = dclus2, scale = TRUE, scores = TRUE) pc <- prcomp(~ api99 + api00 + ell, data = apiclus2, scale = TRUE, retx = TRUE) cov.wt(cbind(svypc$x[, "PC1"], predict(svypc, apiclus2)[, "PC1"]), wt = weights(dclus2), cor = TRUE)$cor # correlação: 0,40 cor(cbind(pc$x[, "PC1"], predict(pc, apiclus2)[, "PC1"])) # correlação: 1 Então, novamente, a divergência entre o $x e o predict() é específica do svyprcomp() e tem a ver com a forma como ele lida com os pesos das observações. Abraços, Leonardo Ferreira Fontenelle[3] Links: 1. http://dx.doi.org/10.5712/rbmfc7%2825%29482 2. https://stat.ethz.ch/pipermail/r-help/2003-July/036645.html 3. http://lattes.cnpq.br/9234772336296638
participantes (2)
-
Cesar Rabak -
Leonardo Ferreira Fontenelle