ANÁLISE FATORIAL

Vamos considerar os exemplos do professor em sala de aula e estudar um pouco sobre o tópico de análise fatorial.

Exemplo 1

  • Considere a matriz de covariância dada abaixo:

\[\Sigma = \begin{pmatrix} 19 & 30 & 2 & 12 \\ 30 & 57 & 5 & 23 \\ 2 & 5 & 38 & 47\\ 12 & 23 & 47 & 68 \end{pmatrix}\]

verifique que \(\Sigma = LL^T + \Psi\) para m = 2, com

\[L = \begin{pmatrix} 4 & 1 \\ 7 & 2 \\ -1 & 6\\ 1 & 8 \end{pmatrix}\]

e \(\Psi = diag(2,4,1,3).\)

SIGMA = matrix(c(19,30,2,12,30,57,5,23,2,5,38,47,12,23,47,68),4,4,byrow = T)
L = matrix(c(4,1,7,2,-1,6,1,8),4,2,byrow = T)
PSI = diag(c(2,4,1,3))
SIGMA
##      [,1] [,2] [,3] [,4]
## [1,]   19   30    2   12
## [2,]   30   57    5   23
## [3,]    2    5   38   47
## [4,]   12   23   47   68
L
##      [,1] [,2]
## [1,]    4    1
## [2,]    7    2
## [3,]   -1    6
## [4,]    1    8
PSI
##      [,1] [,2] [,3] [,4]
## [1,]    2    0    0    0
## [2,]    0    4    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    3

Vamos verificar, se de fato essas cargas fatoriais, com esses dois fatores estão adequados:

L%*%t(L)+PSI # Verificando os valores 
##      [,1] [,2] [,3] [,4]
## [1,]   19   30    2   12
## [2,]   30   57    5   23
## [3,]    2    5   38   47
## [4,]   12   23   47   68

Vamos comparar se é igual a Matriz de covariâncias Sigma.

SIGMA - (L%*%t(L) + PSI)
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0
## [3,]    0    0    0    0
## [4,]    0    0    0    0

Vemos, que de fato isto foi atendido.

  • Verificando as comunalidades.
comunal = function(x) { sum(x^2) }
h = apply(L, 1, "comunal")
h
## [1] 17 53 37 65
h+diag(PSI)
## [1] 19 57 38 68
h+diag(PSI)==diag(SIGMA)
## [1] TRUE TRUE TRUE TRUE

De fato, a soma da comunalidade com a variância específica é igual a variância que se encontra na diagonal da matriz de covariância.

Exemplo 2

  • Para a matriz de covariâncias abaixo, verifique a não existência de solução para m = 1.

\[\Sigma = \begin{pmatrix} 1 & 0.9 & 0.7 \\ 0.9 & 1 & 0.4 \\ 0.7 & 0.4 & 1 \end{pmatrix}\]

# A matriz RHO é positiva definida
RHO    <- matrix(c(1,0.9,0.7,0.9,1,0.4,0.7,0.4,1),3,3,byrow =  T)
RHO
##      [,1] [,2] [,3]
## [1,]  1.0  0.9  0.7
## [2,]  0.9  1.0  0.4
## [3,]  0.7  0.4  1.0

Pode ser obtida uma solução através de componentes principais, todavia ela não é exata, só é exata se tivermos m=p=3, que não haveria necessidade já que o objetivo é diminuir o número de dimensões.

dec    <- eigen(RHO)
delta  <- dec$values #autovalores
P      <- dec$vectors #autovetores
ell.1  <- sqrt(delta[1])*P[,1]
ell.2  <- sqrt(delta[2])*P[,2]
ell.3  <- sqrt(delta[3])*P[,3]
L      <- cbind(ell.1,ell.2,ell.3)
L
##           ell.1      ell.2       ell.3
## [1,] -0.9874197 -0.0871829  0.13191463
## [2,] -0.8846479 -0.4553643 -0.10020701
## [3,] -0.7720339  0.6332923 -0.05389295
RHO-ell.1%*%t(ell.1) # Se considerar m=1 a solucao nao é exata
##             [,1]        [,2]        [,3]
## [1,]  0.02500233  0.02648121 -0.06232153
## [2,]  0.02648121  0.21739805 -0.28297823
## [3,] -0.06232153 -0.28297823  0.40396359
# Aproveitamos o exemplo para verificar:
delta/3         # Proporcao da variabilidade
## [1] 0.78454534 0.20533887 0.01011579
cumsum(delta)/3 # Proporcao da variabilidade acumulada
## [1] 0.7845453 0.9898842 1.0000000
# Para conferir, note que delta_k = ell.1_k'*ell.1_k
# Logo, podemos ver a proporcao da variabilidade devida ao 
# k-esimo fator atraves da soma de quadrados dos 
# valores da k-esima coluna de L
delta
## [1] 2.35363603 0.61601660 0.03034736
cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2,(ell.3)%*%ell.3)
##          [,1]      [,2]       [,3]
## [1,] 2.353636 0.6160166 0.03034736

Exemplo 3

  • Em um estudo sobre as preferências de consumo, solicitou-se a uma amostra de consumidores que classificassem vários atributos de um novo produto. As respostas, em uma escala diferencial semântica de 1 a 7 pontos, foram tabuladas e a matriz de correlações dos atributos foi calculada. Os atributos classificados foram sabor \((X_1)\), preço \((X_2)\), aroma \((X_3)\), “bom como um lanche” \((X_4)\) e “fornece muita energia” \((X_5)\).

\[R = \begin{pmatrix} 1 & 0.02 & 0.96 & 0.42 & 0.01 \\ 0.02 & 1 & 0.13 & 0.71 & 0.85 \\ 0.96 & 0.13 & 1 & 0.50 & 0.11 \\ 0.42 & 0.71 & 0.50 & 1 & 0.79 \\ 0.01 & 0.85 & 0.11 & 0.79 & 1 \end{pmatrix}\]

Pode-se observar que as variáveis 1 e 3 (gosto e aroma) formam um grupo, assim como as variáveis 2 e 5 (preço e energia). A variável 4, “bom para lanche” está mais próxima do grupo (preço, energia).

R <- matrix(c(
  1, 0.02, 0.96, 0.42, 0.01,
  0.02,    1, 0.13, 0.71, 0.85,
  0.96, 0.13,    1, 0.50, 0.11,
  0.42, 0.71, 0.50,    1, 0.79,
  0.01, 0.85, 0.11, 0.79,    1),5,5, byrow = T)
R
##      [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00 0.02 0.96 0.42 0.01
## [2,] 0.02 1.00 0.13 0.71 0.85
## [3,] 0.96 0.13 1.00 0.50 0.11
## [4,] 0.42 0.71 0.50 1.00 0.79
## [5,] 0.01 0.85 0.11 0.79 1.00
  • Realizando a análise fatorial via componentes principais.
p <- 5 # número de variáveis 
AF <- eigen(R)
delta <- AF$values #autovalores
P <- AF$vectors #autovetores
delta
## [1] 2.85309042 1.80633245 0.20449022 0.10240947 0.03367744
round(P,3)
##       [,1]   [,2]   [,3]   [,4]   [,5]
## [1,] 0.331 -0.607  0.098  0.139  0.702
## [2,] 0.460  0.390  0.743 -0.282  0.072
## [3,] 0.382 -0.557  0.168  0.117 -0.709
## [4,] 0.556  0.078 -0.602 -0.568  0.002
## [5,] 0.473  0.404 -0.221  0.751  0.009
delta/p         # Proporcao da variabilidade
## [1] 0.570618084 0.361266491 0.040898045 0.020481893 0.006735487

Verificando a proporção da variabilidade explicada para saber quantos fatores vamos utilizar.

cumsum(delta)/p
## [1] 0.5706181 0.9318846 0.9727826 0.9932645 1.0000000

Como podemos enxergar, com dois fatores temos uma explicação de cerca de 93.19% da variabilidade dos dados. Sendo assim, vamos escolher m = 2.

  • Estimando a matriz de cargas dos fatores, para m = 2.
ell.1 <- sqrt(delta[1])*P[,1]
ell.2 <- sqrt(delta[2])*P[,2]
L     <- cbind(ell.1,ell.2)
L
##          ell.1      ell.2
## [1,] 0.5598618 -0.8160981
## [2,] 0.7772594  0.5242021
## [3,] 0.6453364 -0.7479464
## [4,] 0.9391057  0.1049187
## [5,] 0.7982069  0.5432281
L%*%t(L)
##             [,1]        [,2]      [,3]      [,4]        [,5]
## [1,] 0.979461354 0.007357539 0.9716968 0.4401455 0.003558179
## [2,] 0.007357539 0.878920022 0.1095187 0.7849273 0.905175175
## [3,] 0.971696840 0.109518709 0.9758829 0.5275656 0.108806485
## [4,] 0.440145533 0.784927343 0.5275656 0.8929275 0.806595488
## [5,] 0.003558179 0.905175175 0.1088065 0.8065955 0.932231117
round(R-L%*%t(L),digits = 4)
##         [,1]    [,2]    [,3]    [,4]    [,5]
## [1,]  0.0205  0.0126 -0.0117 -0.0201  0.0064
## [2,]  0.0126  0.1211  0.0205 -0.0749 -0.0552
## [3,] -0.0117  0.0205  0.0241 -0.0276  0.0012
## [4,] -0.0201 -0.0749 -0.0276  0.1071 -0.0166
## [5,]  0.0064 -0.0552  0.0012 -0.0166  0.0678
# Para conferir, note que delta_j = ell.1_j'*ell.1_j
# Logo, podemos ver a proporcao da variabilidade devida ao 
# k-esimo fator atraves da soma de quadrados dos 
# valores da k-esima coluna de L
delta
## [1] 2.85309042 1.80633245 0.20449022 0.10240947 0.03367744
cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2)
##         [,1]     [,2]
## [1,] 2.85309 1.806332
# Estimacao das comunalidades 
sum2 <- function(x){return(sum(x^2))}
h    <- apply(L,1,"sum2")
h
## [1] 0.9794614 0.8789200 0.9758829 0.8929275 0.9322311
# Estimacao das variancias especificas
psi  <- 1 - h # variancia das variaveis padronizadas - comundalidades 
psi
## [1] 0.02053865 0.12107998 0.02411712 0.10707250 0.06776888
R_estimado <- L%*%t(L) + diag(psi)
R_estimado
##             [,1]        [,2]      [,3]      [,4]        [,5]
## [1,] 1.000000000 0.007357539 0.9716968 0.4401455 0.003558179
## [2,] 0.007357539 1.000000000 0.1095187 0.7849273 0.905175175
## [3,] 0.971696840 0.109518709 1.0000000 0.5275656 0.108806485
## [4,] 0.440145533 0.784927343 0.5275656 1.0000000 0.806595488
## [5,] 0.003558179 0.905175175 0.1088065 0.8065955 1.000000000
#matriz residual

U <- R - R_estimado
U
##              [,1]        [,2]         [,3]        [,4]         [,5]
## [1,]  0.000000000  0.01264246 -0.011696840 -0.02014553  0.006441821
## [2,]  0.012642461  0.00000000  0.020481291 -0.07492734 -0.055175175
## [3,] -0.011696840  0.02048129  0.000000000 -0.02756557  0.001193515
## [4,] -0.020145533 -0.07492734 -0.027565574  0.00000000 -0.016595488
## [5,]  0.006441821 -0.05517518  0.001193515 -0.01659549  0.000000000

Obtivemos boa explicação através de dois fatores

  • Aplicando a rotação de fatores, com o método varimax
V.L <- varimax(L)
varimax(V.L$loadings)
## $loadings
## 
## Loadings:
##      ell.1  ell.2 
## [1,]        -0.989
## [2,]  0.937       
## [3,]  0.129 -0.979
## [4,]  0.843 -0.427
## [5,]  0.965       
## 
##                ell.1 ell.2
## SS loadings    2.538 2.121
## Proportion Var 0.508 0.424
## Cumulative Var 0.508 0.932
## 
## $rotmat
##              [,1]          [,2]
## [1,] 0.9999997812 -0.0006614858
## [2,] 0.0006614858  0.9999997812
V.L$loadings
## 
## Loadings:
##      ell.1  ell.2 
## [1,]        -0.989
## [2,]  0.937       
## [3,]  0.130 -0.979
## [4,]  0.843 -0.427
## [5,]  0.965       
## 
##                ell.1 ell.2
## SS loadings    2.539 2.121
## Proportion Var 0.508 0.424
## Cumulative Var 0.508 0.932

Tivemos essas cargas fatoriais e podemos dizer que o primeiro é uma média ponderada, e o segundo também, identificando dois grupos.

cbind(L,V.L$loadings[,1:2])
##          ell.1      ell.2      ell.1       ell.2
## [1,] 0.5598618 -0.8160981 0.02125556 -0.98944912
## [2,] 0.7772594  0.5242021 0.93742129  0.01270189
## [3,] 0.6453364 -0.7479464 0.13009861 -0.97926362
## [4,] 0.9391057  0.1049187 0.84310820 -0.42672715
## [5,] 0.7982069  0.5432281 0.96536898  0.01714214

Exemplo 4

  • Considere, novamente, o exemplo sobre o mercado de ações.
rm(list=ls(all=T))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/ANÁLISE MULTIVARIADA II/AULAS/ANÁLISE FATORIAL")
dados <- read.table("acp-ex4-data.txt", sep = "\t", 
                    dec = ".", header = F,
                    col.names = c("JP", "Citi", "Fargo", "Royal", "Exxon"))
DescTools::Desc(dados, plotit = F, digits = 4)
## ------------------------------------------------------------------------------ 
## Describe dados (data.frame):
## 
## data frame:  103 obs. of  5 variables
##      103 complete cases (100.0%)
## 
##   Nr  ColName  Class    NAs  Levels
##   1   JP       numeric  .          
##   2   Citi     numeric  .          
##   3   Fargo    numeric  .          
##   4   Royal    numeric  .          
##   5   Exxon    numeric  .          
## 
## 
## ------------------------------------------------------------------------------ 
## 1 - JP (numeric)
## 
##    length        n      NAs  unique      0s     mean   meanCI'
##       103      103        0     = n       0   0.0011  -0.0030
##             100.0%     0.0%            0.0%            0.0051
##                                                              
##       .05      .10      .25  median     .75      .90      .95
##   -0.0330  -0.0265  -0.0136  0.0034  0.0168   0.0238   0.0302
##                                                              
##     range       sd    vcoef     mad     IQR     skew     kurt
##    0.0943   0.0208  19.5855  0.0217  0.0304  -0.1128  -0.5367
##                                                              
## lowest : -0.0459, -0.0447, -0.0374, -0.0345, -0.0331
## highest: 0.0308, 0.0373, 0.0466, 0.0485, 0.0485
## 
## ' 95%-CI (classic)
## 
## ------------------------------------------------------------------------------ 
## 2 - Citi (numeric)
## 
##    length        n      NAs  unique      0s    mean   meanCI'
##       103      103        0     = n       0  0.0007  -0.0034
##             100.0%     0.0%            0.0%           0.0047
##                                                             
##       .05      .10      .25  median     .75     .90      .95
##   -0.0311  -0.0227  -0.0132  0.0017  0.0140  0.0239   0.0394
##                                                             
##     range       sd    vcoef     mad     IQR    skew     kurt
##    0.1123   0.0209  31.9575  0.0207  0.0273  0.0893   0.3767
##                                                             
## lowest : -0.0598, -0.0497, -0.0431, -0.0408, -0.0313
## highest: 0.0411, 0.0441, 0.0457, 0.0523, 0.0525
## 
## ' 95%-CI (classic)
## 
## ------------------------------------------------------------------------------ 
## 3 - Fargo (numeric)
## 
##    length        n      NAs  unique      0s    mean   meanCI'
##       103      103        0     102       2  0.0016  -0.0013
##             100.0%     0.0%            1.9%           0.0046
##                                                             
##       .05      .10      .25  median     .75     .90      .95
##   -0.0187  -0.0164  -0.0081  0.0003  0.0100  0.0206   0.0289
##                                                             
##     range       sd    vcoef     mad     IQR    skew     kurt
##    0.0769   0.0150   9.2035  0.0130  0.0181  0.3222   0.0108
##                                                             
## lowest : -0.0362, -0.0304, -0.0238, -0.0238, -0.0220
## highest: 0.0318, 0.0344, 0.0350, 0.0382, 0.0407
## 
## ' 95%-CI (classic)
## 
## ------------------------------------------------------------------------------ 
## 4 - Royal (numeric)
## 
##    length        n      NAs  unique      0s     mean   meanCI'
##       103      103        0     = n       1   0.0040  -0.0012
##             100.0%     0.0%            1.0%            0.0093
##                                                              
##       .05      .10      .25  median     .75      .90      .95
##   -0.0446  -0.0305  -0.0145  0.0063  0.0222   0.0363   0.0445
##                                                              
##     range       sd    vcoef     mad     IQR     skew     kurt
##    0.1159   0.0269   6.6383  0.0245  0.0367  -0.1267  -0.4341
##                                                              
## lowest : -0.0539, -0.0519, -0.0508, -0.0498, -0.0482
## highest: 0.0510, 0.0569, 0.0582, 0.0616, 0.0620
## 
## ' 95%-CI (classic)
## 
## ------------------------------------------------------------------------------ 
## 5 - Exxon (numeric)
## 
##    length        n      NAs  unique      0s     mean   meanCI'
##       103      103        0     = n       0   0.0040  -0.0014
##             100.0%     0.0%            0.0%            0.0094
##                                                              
##       .05      .10      .25  median     .75      .90      .95
##   -0.0445  -0.0264  -0.0125  0.0052  0.0216   0.0381   0.0509
##                                                              
##     range       sd    vcoef     mad     IQR     skew     kurt
##    0.1420   0.0277   6.8515  0.0249  0.0342  -0.1173   0.3017
##                                                              
## lowest : -0.0636, -0.0629, -0.0620, -0.0583, -0.0557
## highest: 0.0561, 0.0574, 0.0587, 0.0649, 0.0784
## 
## ' 95%-CI (classic)
n <- dim(dados)[1] # quantidade de observações
p <- dim(dados)[2] # quantidade de variáveis 
R <- cor(dados) # matriz de correlação das variáveis 
round(R,4)
##           JP   Citi  Fargo  Royal  Exxon
## JP    1.0000 0.6323 0.5105 0.1146 0.1545
## Citi  0.6323 1.0000 0.5741 0.3223 0.2127
## Fargo 0.5105 0.5741 1.0000 0.1825 0.1462
## Royal 0.1146 0.3223 0.1825 1.0000 0.6834
## Exxon 0.1545 0.2127 0.1462 0.6834 1.0000

Análise fatorial via método da máxima verossimilhança

AFMV.1 <- factanal(dados,factors=1,rotation="none") # Com um único fator 
AFMV.2 <- factanal(dados,factors=2,rotation="none") # Com dois fatores

Estimação da matriz de cargas fatoriais

L.1 <- AFMV.1$loadings[,1]
L.2 <- AFMV.2$loadings[,1:2]
L.1
##        JP      Citi     Fargo     Royal     Exxon 
## 0.7154389 0.8750793 0.6630512 0.3460438 0.2792017
L.2
##         Factor1      Factor2
## JP    0.1205972  0.754267060
## Citi  0.3284924  0.785749642
## Fargo 0.1876017  0.650217058
## Royal 0.9974724 -0.007103505
## Exxon 0.6851746  0.026317443

Estimação das comunalidades

sum2  <- function(x){return(sum(x^2))}
h.1   <- L.1^2                 # um fator
h.2   <- apply(L.2,1,"sum2")   # dois fatores
h.1
##         JP       Citi      Fargo      Royal      Exxon 
## 0.51185287 0.76576377 0.43963688 0.11974629 0.07795361
h.2
##        JP      Citi     Fargo     Royal     Exxon 
## 0.5834625 0.7253097 0.4579766 0.9950016 0.4701569

Estimação das variâncias específicas

psi.1    <- AFMV.1$uniqueness
psi.2    <- AFMV.2$uniqueness
psi.1
##        JP      Citi     Fargo     Royal     Exxon 
## 0.4881494 0.2342357 0.5603556 0.8802602 0.9220490
psi.2
##        JP      Citi     Fargo     Royal     Exxon 
## 0.4165374 0.2746902 0.5420233 0.0050000 0.5298429

Estimação/Recomposição de R para m = 1

R.est.1  <- L.1%*%t(L.1)+diag(psi.1)
R.est.1
##             JP      Citi     Fargo      Royal      Exxon
## [1,] 1.0000022 0.6260658 0.4743726 0.24757319 0.19975179
## [2,] 0.6260658 0.9999994 0.5802224 0.30281574 0.24432365
## [3,] 0.4743726 0.5802224 0.9999925 0.22944473 0.18512504
## [4,] 0.2475732 0.3028157 0.2294447 1.00000651 0.09661602
## [5,] 0.1997518 0.2443237 0.1851250 0.09661602 1.00000263

Estimação/Recomposição de R para m = 2

R.est.2  <- L.2%*%t(L.2)+diag(psi.2)
R.est.2
##              JP      Citi     Fargo     Royal     Exxon
## JP    0.9999999 0.6322803 0.5130616 0.1149345 0.1024805
## Citi  0.6322803 1.0000000 0.5725336 0.3220805 0.2457536
## Fargo 0.5130616 0.5725336 0.9999999 0.1825087 0.1456520
## Royal 0.1149345 0.3220805 0.1825087 1.0000016 0.6832558
## Exxon 0.1024805 0.2457536 0.1456520 0.6832558 0.9999997

Matriz Residual para m = 1

U.1 <- R - R.est.1
round(U.1,4)
##            JP    Citi   Fargo   Royal   Exxon
## JP     0.0000  0.0062  0.0361 -0.1330 -0.0453
## Citi   0.0062  0.0000 -0.0061  0.0195 -0.0316
## Fargo  0.0361 -0.0061  0.0000 -0.0469 -0.0389
## Royal -0.1330  0.0195 -0.0469  0.0000  0.5868
## Exxon -0.0453 -0.0316 -0.0389  0.5868  0.0000

Matriz Residual para m = 2

U.2 <- R - R.est.2
round(U.2,4)
##            JP    Citi   Fargo  Royal   Exxon
## JP     0.0000  0.0000 -0.0026 -3e-04  0.0520
## Citi   0.0000  0.0000  0.0016  2e-04 -0.0331
## Fargo -0.0026  0.0016  0.0000  0e+00  0.0006
## Royal -0.0003  0.0002  0.0000  0e+00  0.0001
## Exxon  0.0520 -0.0331  0.0006  1e-04  0.0000

Conforme podemos observar a análise fatorial com dois fatores torna-se mais precisa do que com um, um resultado que já tinha sido observado acima, pois, a proporção da variabilidade explicada por dois fatores foi de: 64.6%.

Soma de quadrados das cargas para o k-ésimo fator

ssc.1 <- sum2(L.1)           # um fator
ssc.2 <- apply(L.2,2,"sum2") # dois fatores
ssc.1
## [1] 1.914953
ssc.2
##  Factor1  Factor2 
## 1.622061 1.609847

Proporção da variabilidade total

ssc.1/p # por um único fator
## [1] 0.3829907
ssc.2/p # por dois fatores
##   Factor1   Factor2 
## 0.3244121 0.3219693
prop.ac <- cumsum(ssc.2)/p # proporção da variabilidade acumulada m = 2
prop.ac 
##   Factor1   Factor2 
## 0.3244121 0.6463815

Teste de adequação

# Testes de adequacao
nu.1 <- ((p-1)^2-p-1)/2
nu.2 <- ((p-2)^2-p-2)/2
nu.1;nu.2
## [1] 5
## [1] 1
m1  <- 1
W1  <- (n-1-(2*p + 4*m1 + 5)/6)*( log(det(R.est.1))-log(det(R)))
W1
## [1] 62.22057

Se o teste abaixo for verificado, podemos continuar que não é adequado a utilização da análise fatorial

W1  > qchisq(0.95,nu.1) 
## [1] TRUE
valorp.1 <- 1-pchisq(W1,nu.1)
valorp.1
## [1] 4.221401e-12

Como a hipótese nula foi rejeitada, a quantidade de um fator não está adequada

m2  <- 2
W2  <- (n-1-(2*p + 4* m2 + 5)/6)*( log(det(R.est.2))-log(det(R)))
W2
## [1] 2.004672
W2  > qchisq(0.95,nu.2) 
## [1] FALSE
valorp.2 <- 1-pchisq(W2,nu.2)
valorp.2
## [1] 0.1568152

Como a hipótese nula não foi rejeitada, a quantidade de dois fatores está adequada.

AFMV.2$PVAL
## objective 
## 0.1600081

Comparado com o teste, vemos que dois fatores está bem adequado. Com um único fator, o resultado do p-valor foi de 4.22^{-12}. #### Estimação dos Escores dos fatores

AFMV.2.Reg <- factanal(dados,factors=2,rotation="none",scores="regression")
AFMV.2.Bar <- factanal(dados,factors=2,rotation="none",scores="Bartlett")
  • Plotando os gráficos
plot(AFMV.2.Reg$scores[,1],AFMV.2.Bar$scores[,2])

plot(AFMV.2.Reg$scores[,2],AFMV.2.Bar$scores[,1])

cor(AFMV.2.Reg$scores, AFMV.2.Bar$scores)
##              Factor1      Factor2
## Factor1 1.000000e+00 1.840696e-17
## Factor2 3.429087e-17 1.000000e+00

Análise fatorial via componentes principais

AF      <- prcomp(dados,scale=T)
delta   <- AF$sdev^2
P       <- AF$rotation 
delta
## [1] 2.4372731 1.4070127 0.5005127 0.4000316 0.2551699
P
##              PC1        PC2         PC3        PC4         PC5
## JP    -0.4690832  0.3680070 -0.60431522  0.3630228  0.38412160
## Citi  -0.5324055  0.2364624 -0.13610618 -0.6292079 -0.49618794
## Fargo -0.4651633  0.3151795  0.77182810  0.2889658  0.07116948
## Royal -0.3873459 -0.5850373  0.09336192 -0.3812515  0.59466408
## Exxon -0.3606821 -0.6058463 -0.10882629  0.4934145 -0.49755167

Proporção da variabilidade para m = 2, dois componentes principais

delta/p
## [1] 0.48745462 0.28140253 0.10010255 0.08000632 0.05103398
cumsum(delta)/p
## [1] 0.4874546 0.7688572 0.8689597 0.9489660 1.0000000

Ou seja, com dois fatores através de componentes principais temos: 76.89%.

Estimação da matriz de cargas fatoriais

ell.1 <- sqrt(delta[1])*P[,1]
ell.2 <- sqrt(delta[2])*P[,2]
L.1   <- ell.1               # um fator
L.2   <- cbind(ell.1,ell.2)  # dois fatores
L.1
##         JP       Citi      Fargo      Royal      Exxon 
## -0.7323218 -0.8311791 -0.7262022 -0.6047155 -0.5630885
L.2
##            ell.1      ell.2
## JP    -0.7323218  0.4365209
## Citi  -0.8311791  0.2804859
## Fargo -0.7262022  0.3738582
## Royal -0.6047155 -0.6939569
## Exxon -0.5630885 -0.7186401

Conforme já tinha sido encontrado através dos componentes principais estes dois escores, o primeiro uma soma ponderada e o segundo um contraste entre as ações dos bancos e as de petróleo.

Estimação das comunalidades

h.1   <- L.1^2                 # um fator
h.2   <- apply(L.2,1,"sum2")   # dois fatores
h.1
##        JP      Citi     Fargo     Royal     Exxon 
## 0.5362953 0.6908587 0.5273697 0.3656808 0.3170686
h.2
##        JP      Citi     Fargo     Royal     Exxon 
## 0.7268458 0.7695311 0.6671396 0.8472571 0.8335122

Estimação das variânciais específicas

psi_1 <- 1 - h.1
psi_2 <- 1 - h.2
psi_1
##        JP      Citi     Fargo     Royal     Exxon 
## 0.4637047 0.3091413 0.4726303 0.6343192 0.6829314
psi_2
##        JP      Citi     Fargo     Royal     Exxon 
## 0.2731542 0.2304689 0.3328604 0.1527429 0.1664878

Estimação/Recomposição de R para m = 2

R.est_1  <- L.1%*%t(L.1)+diag(psi_1)
R.est_2  <- L.2%*%t(L.2)+diag(psi_2)
R.est_1
##             JP      Citi     Fargo     Royal     Exxon
## [1,] 1.0000000 0.6086906 0.5318137 0.4428464 0.4123620
## [2,] 0.6086906 1.0000000 0.6036041 0.5026269 0.4680274
## [3,] 0.5318137 0.6036041 1.0000000 0.4391457 0.4089161
## [4,] 0.4428464 0.5026269 0.4391457 1.0000000 0.3405083
## [5,] 0.4123620 0.4680274 0.4089161 0.3405083 1.0000000
R.est_2
##               JP      Citi     Fargo     Royal      Exxon
## JP    1.00000000 0.7311286 0.6950107 0.1399197 0.09866056
## Citi  0.73112857 1.0000000 0.7084661 0.3079818 0.26645897
## Fargo 0.69501067 0.7084661 1.0000000 0.1797042 0.14024657
## Royal 0.13991966 0.3079818 0.1797042 1.0000000 0.83921362
## Exxon 0.09866056 0.2664590 0.1402466 0.8392136 1.00000000

Matriz Residual

  • Onde podemos ver se ficou bom a análise fatorial
U_1 <- R - R.est_2
round(U_1,3)
##           JP   Citi  Fargo  Royal  Exxon
## JP     0.000 -0.099 -0.185 -0.025  0.056
## Citi  -0.099  0.000 -0.134  0.014 -0.054
## Fargo -0.185 -0.134  0.000  0.003  0.006
## Royal -0.025  0.014  0.003  0.000 -0.156
## Exxon  0.056 -0.054  0.006 -0.156  0.000