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

Terça Abril 26 23:56:07 BRT 2016

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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160426/84b5bf98/attachment.html>