#’ libraries needed

library(rio)
library(magrittr)
library(polycor)
library(psych)
## 
## Adjuntando el paquete: 'psych'
## The following object is masked from 'package:polycor':
## 
##     polyserial
library(matrixcalc)
library(GPArotation)
## 
## Adjuntando el paquete: 'GPArotation'
## The following objects are masked from 'package:psych':
## 
##     equamax, varimin
library(BBmisc)
## 
## Adjuntando el paquete: 'BBmisc'
## The following object is masked from 'package:base':
## 
##     isFALSE
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:BBmisc':
## 
##     coalesce, collapse, symdiff
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#’ read data

data_preg1=import("reporte_com.xlsx")
names(data_preg1)
## [1] "PROVINCIA"  "COM_NO_ELE" "COM_SI_ELE" "COM_NO_GLP" "COM_SI_GLP"
## [6] "COM_NO_CAR" "COM_SI_CAR" "COM_NO_LE"  "COM_SI_LE"
str(data_preg1)
## 'data.frame':    196 obs. of  9 variables:
##  $ PROVINCIA : chr  "ABANCAY" "ACOBAMBA" "ACOMAYO" "AIJA" ...
##  $ COM_NO_ELE: num  33780 11600 7037 1989 28563 ...
##  $ COM_SI_ELE: num  490 70 118 16 446 152 447 87 127 31 ...
##  $ COM_NO_GLP: num  10121 8253 4398 1527 15387 ...
##  $ COM_SI_GLP: num  24149 3417 2757 478 13622 ...
##  $ COM_NO_CAR: num  34111 11623 7113 2001 26737 ...
##  $ COM_SI_CAR: num  159 47 42 4 2272 ...
##  $ COM_NO_LE : num  16860 1878 841 219 11404 ...
##  $ COM_SI_LE : num  17410 9792 6314 1786 17605 ...

#’ data for factorial

dataporc <- data_preg1 %>%
  mutate(
    PORC_SI_ELE = (COM_SI_ELE / (COM_SI_ELE + COM_NO_ELE)) * 100,  # Porcentaje de "Sí electricidad"
    PORC_NO_ELE = (COM_NO_ELE / (COM_SI_ELE + COM_NO_ELE)) * 100,   # Porcentaje de "No electricidad"
    
    PORC_SI_GLP = (COM_SI_GLP / (COM_SI_GLP + COM_NO_GLP)) * 100,  # Porcentaje de "Sí glp"
    PORC_NO_GLP = (COM_NO_GLP / (COM_SI_GLP + COM_NO_GLP)) * 100,   # Porcentaje de "No glp"
    
    PORC_SI_CAR = (COM_SI_CAR / (COM_SI_CAR + COM_NO_CAR)) * 100,  
    PORC_NO_CAR = (COM_NO_CAR / (COM_SI_CAR + COM_NO_CAR)) * 100,
    
    PORC_SI_LE = (COM_SI_LE / (COM_SI_LE + COM_NO_LE)) * 100,  
    PORC_NO_LE = (COM_NO_LE / (COM_SI_LE + COM_NO_LE)) * 100
    
  ) %>%
  select(PROVINCIA, PORC_SI_ELE, PORC_NO_ELE, PORC_SI_GLP, PORC_NO_GLP, PORC_SI_CAR, PORC_NO_CAR, PORC_SI_LE, PORC_NO_LE) 

#’ select data

dontselect=c("PROVINCIA", "PORC_NO_ELE", "PORC_NO_GLP","PORC_NO_CAR","PORC_NO_LE")
select=setdiff(names(dataporc),dontselect) 
thedata_preg1=dataporc[,select]

head(thedata_preg1,10)%>%
  rmarkdown::paged_table()

#’ correlations

corMatrix=polycor::hetcor(thedata_preg1)$correlations

#’ previous evaluations

KMO(corMatrix)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = corMatrix)
## Overall MSA =  0.63
## MSA for each item = 
## PORC_SI_ELE PORC_SI_GLP PORC_SI_CAR  PORC_SI_LE 
##        0.91        0.59        0.72        0.59
cortest.bartlett(corMatrix,n=nrow(thedata_preg1))$p.value>0.05
## [1] FALSE
is.singular.matrix(corMatrix)
## [1] FALSE

#’ number of factors

fa.parallel(thedata_preg1, fa = 'fa',correct = T,plot = F)
## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

#’ varimax

resfa1 <- fa(thedata_preg1,
               nfactors = 1,
               cor = 'mixed',
               rotate = "varimax", #oblimin?
               fm="minres")
print(resfa1$loadings)
## 
## Loadings:
##             MR1   
## PORC_SI_ELE  0.491
## PORC_SI_GLP  0.912
## PORC_SI_CAR  0.280
## PORC_SI_LE  -0.956
## 
##                  MR1
## SS loadings    2.065
## Proportion Var 0.516

#’ oblimin

resfa2 <- fa(thedata_preg1,
               nfactors = 1,
               cor = 'mixed',
               rotate = "oblimin", #oblimin?
               fm="minres")
print(resfa2$loadings)
## 
## Loadings:
##             MR1   
## PORC_SI_ELE  0.491
## PORC_SI_GLP  0.912
## PORC_SI_CAR  0.280
## PORC_SI_LE  -0.956
## 
##                  MR1
## SS loadings    2.065
## Proportion Var 0.516