IMPORTAR DATAS

library(rio)
datita = import("reporte.xlsx")

## New names:
## • `` -> `...1`
## • `Total` -> `Total...6`
## • `Total` -> `Total...9`
## • `Total` -> `Total...12`
## • `Total` -> `Total...15`

str(datita)

## 'data.frame':    199 obs. of  15 variables:
##  $ ...1                  : chr  NA NA NA NA ...
##  $ Código                : num  101 102 103 104 105 106 107 201 202 203 ...
##  $ Provincia             : chr  "Amazonas, provincia: Chachapoyas" "Amazonas, provincia: Bagua" "Amazonas, provincia: Bongara" "Amazonas, provincia: Condorcanqui" ...
##  $ No usa electricidad   : num  14763 20313 7689 9853 13112 ...
##  $ Sí usa electricidad   : num  574 161 124 14 90 65 255 921 16 33 ...
##  $ Total...6             : num  15337 20474 7813 9867 13202 ...
##  $ No usa gas (balón GLP): num  4696 10557 3154 8331 6863 ...
##  $ Sí usa gas (balón GLP): num  10641 9917 4659 1536 6339 ...
##  $ Total...9             : num  15337 20474 7813 9867 13202 ...
##  $ No usa carbón         : num  15161 20185 7755 9841 13169 ...
##  $ Sí usa carbón         : num  176 289 58 26 33 26 335 218 4 4 ...
##  $ Total...12            : num  15337 20474 7813 9867 13202 ...
##  $ No usa leña           : num  7236 7357 2345 1059 1833 ...
##  $ Sí usa leña           : num  8101 13117 5468 8808 11369 ...
##  $ Total...15            : num  15337 20474 7813 9867 13202 ...

names(datita)

##  [1] "...1"                   "Código"                 "Provincia"             
##  [4] "No usa electricidad"    "Sí usa electricidad"    "Total...6"             
##  [7] "No usa gas (balón GLP)" "Sí usa gas (balón GLP)" "Total...9"             
## [10] "No usa carbón"          "Sí usa carbón"          "Total...12"            
## [13] "No usa leña"            "Sí usa leña"            "Total...15"

ELECTRICIDAD

colnames(datita)[colnames(datita) == "No usa electricidad"] <- "elec_no"

colnames(datita)[colnames(datita) == "Sí usa electricidad"] <- "elec_si"

colnames(datita)[colnames(datita) == "Total...6"] <- "elec_total"

datita$elec_percent <- (datita$elec_si / datita$elec_total) * 100

GAS

colnames(datita)[colnames(datita) == "No usa gas (balón GLP)"] <- "gas_no"

colnames(datita)[colnames(datita) == "Sí usa gas (balón GLP)"] <- "gas_si"

colnames(datita)[colnames(datita) == "Total...9"] <- "gas_total"

datita$gas_percent <- (datita$gas_si / datita$gas_total) * 100

CARBON

colnames(datita)[colnames(datita) == "No usa carbón"] <- "carbon_no"

colnames(datita)[colnames(datita) == "Sí usa carbón"] <- "carbon_si"

colnames(datita)[colnames(datita) == "Total...12"] <- "carbon_total"

datita$carbon_percent <- (datita$carbon_si / datita$carbon_total) * 100

LEÑA

colnames(datita)[colnames(datita) == "No usa leña"] <- "leña_no"

colnames(datita)[colnames(datita) == "Sí usa leña"] <- "leña_si"

colnames(datita)[colnames(datita) == "Total...15"] <- "leña_total"

datita$leña_percent <- (datita$leña_si / datita$leña_total) * 100

ANÁLISIS FACTORIAL

dontselect=c("...1","Código","Provincia",
             "elec_no","elec_si","elec_total",
             "gas_si", "gas_no", "gas_total",
             "leña_si", "leña_no", "leña_total",
             "carbon_si", "carbon_no", "carbon_total")
select=setdiff(names(datita),dontselect) 
theData=datita[,select]

# usaremos:
library(magrittr)
head(theData,10)%>%
    rmarkdown::paged_table()

library(polycor)

corMatrix=polycor::hetcor(theData)$correlations

library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:polycor':
## 
##     polyserial

psych::KMO(corMatrix) 

## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = corMatrix)
## Overall MSA =  0.63
## MSA for each item = 
##   elec_percent    gas_percent carbon_percent   leña_percent 
##           0.91           0.59           0.72           0.59

cortest.bartlett(corMatrix,n=nrow(theData))$p.value>0.05

## [1] FALSE

library(matrixcalc)

is.singular.matrix(corMatrix)

## [1] FALSE

round(corMatrix,2)

##                elec_percent gas_percent carbon_percent leña_percent
## elec_percent           1.00        0.47           0.10        -0.46
## gas_percent            0.47        1.00           0.23        -0.87
## carbon_percent         0.10        0.23           1.00        -0.31
## leña_percent          -0.46       -0.87          -0.31         1.00

fa.parallel(theData, fa = 'fa',correct = T,plot = F)

## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

library(GPArotation)

## 
## Attaching package: 'GPArotation'

## The following objects are masked from 'package:psych':
## 
##     equamax, varimin

resfa <- fa(theData,
            nfactors = 1,
            cor = 'mixed',
            rotate = "varimax", #oblimin?
            fm="minres")
print(resfa$loadings)

## 
## Loadings:
##                MR1   
## elec_percent    0.491
## gas_percent     0.912
## carbon_percent  0.280
## leña_percent   -0.956
## 
##                  MR1
## SS loadings    2.065
## Proportion Var 0.516

library(GPArotation)
resfa2 <- fa(theData,
            nfactors = 1,
            cor = 'mixed',
            rotate = "oblimin", #oblimin?
            fm="minres")
print(resfa2$loadings)

## 
## Loadings:
##                MR1   
## elec_percent    0.491
## gas_percent     0.912
## carbon_percent  0.280
## leña_percent   -0.956
## 
##                  MR1
## SS loadings    2.065
## Proportion Var 0.516

sort(resfa$communality)

## carbon_percent   elec_percent    gas_percent   leña_percent 
##     0.07838591     0.24104291     0.83166545     0.91346690

resfa$TLI

## [1] 0.9605126

resfa$rms

## [1] 0.03000416

resfa$BIC

## [1] -4.150711

resfa$RMSEA

##      RMSEA      lower      upper confidence 
## 0.10545186 0.01954305 0.20220497 0.90000000