# Abrir la matriz de datos RMB1986.rda
load(url("http://pagines.uab.cat/plopez/sites/pagines.uab.cat.plopez/files/RMB1986.rda"))

# install.packages("Hmisc")
library(Hmisc) # Para describe, rcorr
# install.packages("Deducer")
library(Deducer) # Para descriptive.table
# install.packages("corrplot")
library(corrplot) # Para corrplot
# install.packages("psych")
library(psych) # Para cor.plot, cortest.bartlett, principal

# Selección de las variables
attach(RMB1986)
DF=data.frame(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,
              parados_ocupados,parados_1erempleo,mujeresactivas,profesion_alta,profesion_baja,
              terciario_medio, terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio)

# Descriptivos
descriptive.table(vars=d(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,
                         parados_ocupados,parados_1erempleo,mujeresactivas,profesion_alta,profesion_baja,
                         terciario_medio, terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio),
                  data=DF,func.names=c("Mean","St. Deviation","Valid N"))
describe(DF)

########  1. Análisis de la matriz de correlaciones

# Matriz de correlaciones, determinante y gráfico
R=cor(DF,use="pairwise.complete.obs")
R
round(R,digits=1)
det(R)

# Significación de las correlaciones
rcorr(as.matrix(DF), type="pearson")
R2=cor.matrix(variables=d(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,
                         parados_ocupados,parados_1erempleo,mujeresactivas,profesion_alta,profesion_baja,
                         terciario_medio, terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio), 
             data=RMB1986, test=cor.test, method='pearson', alternative="two.sided")
print(R2,N=FALSE,CI=FALSE,stat=FALSE)

# Gráficos de correlaciones
cor.plot(R)
corrplot(R, method="circle")
# qscatter_array(d(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,parados_ocupados,parados_1erempleo,
#                  mujeresactivas,profesion_alta,profesion_baja,terciario_medio,terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio),
#                d(menores,mayores,inmigracion,llegados,analfabetos,univesitarios,escolarizacion14_24,parados_ocupados,parados_1erempleo,
#                  mujeresactivas,profesion_alta,profesion_baja,terciario_medio,terciario_alto,agropecuario,vehiculotrabajo,seccion_municipio),
#                data=RMB1986) + geom_smooth(method="lm") # Diagramas de dispersión en forma matricial

# KMO y Medida de adecuación muestral (se precisa la librería psych)
KMO(R)
# Test de Bartlett. Es necesario especificar el nº de casos de la matriz: "n="
cortest.bartlett(R,n=3509)

########  2. ACP sin rotación

# Análisis de componentes con el paquete psych
ACP1=principal(DF, nfactors= 4, rotate="none", covar=FALSE, scores=TRUE)
ACP1
# Ordenación de la matriz de componentes o factores de carga
print(ACP1,sort=T)
fa.diagram(ACP1)
fa.parallel(DF)

# Gráficos factoriales
biplot(ACP1)
biplot(ACP1,choose=c(1,2))
biplot(ACP1,choose=c(3,4))
etiquetas=c("menores","mayores","inmigracion","llegados","analfabetos","univesitarios",
            "escolarizacion14_24","parados_ocupados","parados_1erempleo","mujeresactivas",
            "profesion_alta","profesion_baja","terciario_medio","terciario_alto","agropecuario",
            "vehiculotrabajo","seccion_municipio")
factor.plot(ACP1,labels=etiquetas,choose=c(1,2))
factor.plot(ACP1,labels=etiquetas,choose=c(3,4))

########  3. ACP con rotación varimax

ACP2=principal(DF, nfactors= 4, rotate="varimax", covar=FALSE, scores=TRUE)
print(ACP2,sort=T)
fa.diagram(ACP2)
fa.parallel(DF)
# Gráficos factoriales
biplot(ACP2)
biplot(ACP2,choose=c(1,2))
biplot(ACP2,choose=c(3,4))
factor.plot(ACP2,labels=etiquetas,choose=c(1,2))
factor.plot(ACP2,labels=etiquetas,choose=c(3,4))

# Para guardar las variables factoriales junto a las originales
RMB1986F1=data.frame(RMB1986,ACP1$scores) # sin rotación
save(RMB1986F1,file='C:/Datos/RMB1986F1.rda') # Es preciso cambiar la carpeta donde se guarda
RMB1986F2=data.frame(RMB1986,ACP2$scores) # con rotación
save(RMB1986F2,file='C:/Datos/RMB1986F2.rda') # Es preciso cambiar la carpeta donde se guarda


#########################################

# Análisis de componentes con el paquete FactoMineR
# install.packages("FactoMineR")
library(FactoMineR)
ACP=PCA(DF, ncp=4, graph=F)
ACP

# Descriptivos
ACP$call$centre # media
ACP$call$ecart.type # desviación

# Valores propios y varianza explicada
ACP$eig

# Coordenadas, correlaciones y contribuciones entre las variables originales y los factores
ACP$var
# Las anteriores por separado y redondeando a 3 decimales
round(ACP$var$coord,digits=3)
round(ACP$var$cos2,digits=3)
round(ACP$var$contrib,digits=3)

# Correlación entre variables y factores, significación
dimdesc(ACP, axes=c(1,2))

# Coordenadas, correlaciones y contribuciones entre los individuos y los factores
head(ACP$ind$coord)
head(ACP$ind$cos2)
head(ACP$ind$contrib)

# Gráficos factoriales
plot.PCA(ACP, axes=c(1,2), choix="var")
plot.PCA(ACP, axes=c(3,4), choix="var")
plot.PCA(ACP, axes=c(1,2), choix="ind")
plot.PCA(ACP, axes=c(3,4), choix="ind")
RMB1986F=data.frame(RMB1986,ACP$ind)

# Rotación varimax, stats package
varimax(ACP$var$coord, normalize=T)




