In [1]:
# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Correlacao amostral
#==============================================================================
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,T)
In [2]:
#==============================================================================
# Analise fatorial via componentes principais
#==============================================================================
p       <- 5 # numero de variaveis
AF      <- eigen(R)
delta   <- AF$values
P       <- AF$vector
print("==== Autovalores ====")
print(delta)
print("==== Autovetores ====")
print(P)
print("==== Proporcao da variabilidade ====")
print(delta/p)
print("==== Proporcao da variabilidade acumulada ====")
print(cumsum(delta)/p)
# Escolhemos m=2

[1] "==== Autovalores ===="
[1] 2,85309042 1,80633245 0,20449022 0,10240947 0,03367744
[1] "==== Autovetores ===="
          [,1]        [,2]        [,3]       [,4]         [,5]
[1,] 0,3314539 -0,60721643 -0,09848524  0,1386643  0,701783012
[2,] 0,4601593  0,39003172 -0,74256408 -0,2821170  0,071674637
[3,] 0,3820572 -0,55650828 -0,16840896  0,1170037 -0,708716714
[4,] 0,5559769  0,07806457  0,60158211 -0,5682357  0,001656352
[5,] 0,4725608  0,40418799  0,22053713  0,7513990  0,009012569
[1] "==== Proporcao da variabilidade ===="
[1] 0,570618084 0,361266491 0,040898045 0,020481893 0,006735487
[1] "==== Proporcao da variabilidade acumulada ===="
[1] 0,5706181 0,9318846 0,9727826 0,9932645 1,0000000
In [3]:
#==============================================================================
# Estimacao da 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)
print(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
In [4]:
#==============================================================================
# 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
#==============================================================================
print(delta)
print(cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2))

[1] 2,85309042 1,80633245 0,20449022 0,10240947 0,03367744
        [,1]     [,2]
[1,] 2,85309 1,806332
In [5]:
#==============================================================================
# Estimacao das comunalidades 
#==============================================================================
sum2 <- function(x){return(sum(x^2))}
h    <- apply(L,1,"sum2")
print(h)

[1] 0,9794614 0,8789200 0,9758829 0,8929275 0,9322311
In [6]:
#==============================================================================
# Estimacao das variancias especificas
#==============================================================================
psi  <- 1 - h # variancia das variaveis padronizadas - comundalidades 
print(psi)

[1] 0,02053865 0,12107998 0,02411712 0,10707250 0,06776888
In [7]:
#==============================================================================
# Estimacao/Recomposicao de R
#==============================================================================
R.est  <- L%*%t(L)+diag(psi)
print(R.est)

            [,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
In [8]:
#==============================================================================
# Matriz residual
#==============================================================================
U <- R - R.est
print(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
In [9]:
#==============================================================================
# Rotacao de fatores
#==============================================================================
print(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
In [10]:
#==============================================================================
# Utilizaremos o metodo varimax
#==============================================================================
V.L <- varimax(L)
print(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

$rotmat
          [,1]       [,2]
[1,] 0,8365697 -0,5478605
[2,] 0,5478605  0,8365697
In [11]:
#==============================================================================
# Grafico da rotacao dos fatores
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.3,cex.axis=1.3,
     lab=c(7,7,5),mar=c(4.5,5.5,4,2.5),xpd=F,cex.main=2.0,bty="n")
plot(L,pch=15,xlab=expression(hat(f)[1]),ylab=expression(hat(f)[2]),
    xlim=c(-1,0),ylim=c(-0.6,0.8))
lines(c(-V.L$rotmat[1,1],0),c(-V.L$rotmat[2,1],0),col="red",lwd=2)
lines(c( V.L$rotmat[1,2],0),c( V.L$rotmat[2,2],0),col="blue",lwd=2)
abline(v=0,col="gray",lwd=2,lty=3)
abline(h=0,col="gray",lwd=2,lty=3)
text(-0.58, 0.78,"1")
text(-0.75,-0.54,"2")
text(-0.68, 0.75,"3")
text(-0.91,-0.10,"4")
text(-0.79,-0.59,"5")
text(-0.45,-0.20,expression(hat(f)[1]^"*"),cex=1.3)
text(-0.45, 0.60,expression(hat(f)[2]^"*"),cex=1.3)
points(-L,pch=15)
No description has been provided for this image
In [12]:
#==============================================================================
# Poderiamos trocar os sinais correspondente ao primeiro fator
# Utilizariamos outra rotacao
#==============================================================================
Q   <- matrix(c(-1,0,0,1),2,2,T)
print(Q)
print(V.L$loadings[,1:2]%*%Q)

     [,1] [,2]
[1,]   -1    0
[2,]    0    1
            [,1]        [,2]
[1,] -0,02125556 -0,98944912
[2,] -0,93742129  0,01270189
[3,] -0,13009861 -0,97926362
[4,] -0,84310820 -0,42672715
[5,] -0,96536898  0,01714214