Entregable 3 - Gabriela Eliana Rollano Málaga
2022-06-20
library(rio)
library(DescTools)
library(ggplot2)
library(moments)
library(Rmisc)## Loading required package: lattice## Loading required package: plyrlibrary(e1071)##
## Attaching package: 'e1071'## The following objects are masked from 'package:moments':
##
## kurtosis, moment, skewnesslibrary(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following objects are masked from 'package:DescTools':
##
## AUC, ICC, SDlibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(gplots)## Registered S3 method overwritten by 'gplots':
## method from
## reorder.factor DescTools##
## Attaching package: 'gplots'## The following object is masked from 'package:DescTools':
##
## reorder.factor## The following object is masked from 'package:stats':
##
## lowesslibrary(vcd)## Loading required package: gridlibrary(PMCMRplus)
library(nortest)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logit## The following object is masked from 'package:DescTools':
##
## Recodelibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(lm.beta)
library(gtools)##
## Attaching package: 'gtools'## The following object is masked from 'package:car':
##
## logit## The following object is masked from 'package:psych':
##
## logit## The following object is masked from 'package:e1071':
##
## permutationslibrary(jtools)##
## Attaching package: 'jtools'## The following object is masked from 'package:DescTools':
##
## %nin%library(ggstance)##
## Attaching package: 'ggstance'## The following objects are masked from 'package:ggplot2':
##
## geom_errorbarh, GeomErrorbarhlibrary(fastDummies)
library(writexl)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(polycor)##
## Attaching package: 'polycor'## The following object is masked from 'package:psych':
##
## polyseriallibrary(ggcorrplot)
library(matrixcalc)
library(GPArotation)
library(lavaan)## This is lavaan 0.6-11
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covlibrary(BBmisc)##
## Attaching package: 'BBmisc'## The following object is masked from 'package:jtools':
##
## %nin%## The following object is masked from 'package:grid':
##
## explode## The following objects are masked from 'package:dplyr':
##
## coalesce, collapse## The following object is masked from 'package:DescTools':
##
## %nin%## The following object is masked from 'package:base':
##
## isFALSEvdem = import("https://github.com/GabrielaRollano/Entregable-3/blob/main/V-Dem-CY-Core-v12.rds?raw=true")
View(vdem)#1. Variable dependiente - Índice de Democracia Liberal (v2x_libdem)
summary(vdem$v2x_libdem)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0050 0.0580 0.1250 0.2209 0.2920 0.8960 2619str(vdem$v2x_libdem)## num [1:27380] 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 ...#2. Variables independientes - Gabriela Rollano #2.1.Autonomía de los Organismos Electorales (v2elembaut_ord)
str(vdem$v2elembaut_ord)## num [1:27380] 0 0 0 0 0 0 0 0 0 0 ...summary(vdem$v2elembaut_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 0.000 1.227 3.000 4.000 3016#2.2 Medios de comunicación escritos o transmitidos críticos del gobierno (v2mecrit_ord)
str(vdem$v2mecrit_ord)## num [1:27380] 1 1 1 1 1 1 1 1 1 1 ...summary(vdem$v2mecrit_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 1.000 1.000 1.397 2.000 3.000 873#2.3. Respeto a Contraargumentos (v2dlcountr_ord)
str(vdem$v2dlcountr_ord)## num [1:27380] NA NA NA NA NA NA NA NA NA NA ...summary(vdem$v2dlcountr_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 1.0 2.0 1.6 2.0 5.0 8370#3. Variables independientes #3.1. Organizaciones de Sociedad civil #3.1.1. Represión gubernamental de Organizaciones de Sociedad Civil (v2csreprss_ord)
str(vdem$v2csreprss_ord)## num [1:27380] 1 1 1 1 1 1 1 1 1 1 ...summary(vdem$v2csreprss_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 1.00 2.00 2.07 3.00 4.00 749#3.1.2. Control gubernamental sobre Organizaciones de Sociedad Civil (v2cseeorgs_ord)
str(vdem$v2cseeorgs_ord)## num [1:27380] 0 0 0 0 0 0 0 0 0 0 ...summary(vdem$v2cseeorgs_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 1.000 1.597 3.000 4.000 580#3.1.3. Consulta gubernamental hacia Organizaciones de Sociedad Civil (v2cscnsult_ord)
str(vdem$v2cscnsult_ord)## num [1:27380] 0 0 0 0 0 0 0 0 0 0 ...summary(vdem$v2cscnsult_ord)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 0.0000 0.4807 1.0000 2.0000 554I. CLUSTERIZACIÓN
Armar la base de apoyo
cluster_gabriela = subset(vdem, select = c(country_name, year, v2x_libdem, v2elembaut_ord, v2mecrit_ord, v2dlcountr_ord, v2csreprss_ord, v2cseeorgs_ord, v2cscnsult_ord))cluster_gabriela = cluster_gabriela[cluster_gabriela$year==2021,]# Trabajaremos con los nombres de los países y las seis variables independientes escogidas, por eso anulamos el año y la variable dependiente.
cluster_gabriela$year = NULL
cluster_gabriela$v2x_libdem = NULLstr(cluster_gabriela)## 'data.frame': 179 obs. of 7 variables:
## $ country_name : chr "Mexico" "Suriname" "Sweden" "Switzerland" ...
## $ v2elembaut_ord: num 4 3 4 4 3 3 3 0 0 2 ...
## $ v2mecrit_ord : num 2 3 3 3 2 3 3 0 1 2 ...
## $ v2dlcountr_ord: num 2 2 4 4 2 3 4 1 1 2 ...
## $ v2csreprss_ord: num 3 4 4 4 4 3 4 1 1 4 ...
## $ v2cseeorgs_ord: num 3 4 4 4 4 4 4 0 1 3 ...
## $ v2cscnsult_ord: num 0 1 2 2 2 1 1 0 0 1 ...Hacemos la estandarización.
library(BBmisc)
cluster_gabriela[,-1]=normalize(cluster_gabriela[,-1],method='standardize')
cluster_gabriela=cluster_gabriela[complete.cases(cluster_gabriela),]summary(cluster_gabriela)## country_name v2elembaut_ord v2mecrit_ord v2dlcountr_ord
## Length:179 Min. :-1.5482 Min. :-2.34382 Min. :-1.9144
## Class :character 1st Qu.:-0.8571 1st Qu.:-0.04495 1st Qu.:-0.1156
## Mode :character Median : 0.5251 Median :-0.04495 Median :-0.1156
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.5251 3rd Qu.: 1.10448 3rd Qu.: 0.7838
## Max. : 1.2162 Max. : 1.10448 Max. : 2.5826
## v2csreprss_ord v2cseeorgs_ord v2cscnsult_ord
## Min. :-2.2807 Min. :-2.0256 Min. :-1.45913
## 1st Qu.:-0.6445 1st Qu.:-0.8459 1st Qu.: 0.03335
## Median : 0.1737 Median : 0.3339 Median : 0.03335
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.9918 3rd Qu.: 1.1205 3rd Qu.: 0.03335
## Max. : 0.9918 Max. : 1.1205 Max. : 1.52583Vemos las correlaciones.
cor(cluster_gabriela[,-1])## v2elembaut_ord v2mecrit_ord v2dlcountr_ord v2csreprss_ord
## v2elembaut_ord 1.0000000 0.7377899 0.6826106 0.7945379
## v2mecrit_ord 0.7377899 1.0000000 0.6917230 0.7422000
## v2dlcountr_ord 0.6826106 0.6917230 1.0000000 0.7146772
## v2csreprss_ord 0.7945379 0.7422000 0.7146772 1.0000000
## v2cseeorgs_ord 0.7581324 0.7667901 0.6945539 0.8888291
## v2cscnsult_ord 0.6429916 0.6375932 0.6675044 0.7144543
## v2cseeorgs_ord v2cscnsult_ord
## v2elembaut_ord 0.7581324 0.6429916
## v2mecrit_ord 0.7667901 0.6375932
## v2dlcountr_ord 0.6945539 0.6675044
## v2csreprss_ord 0.8888291 0.7144543
## v2cseeorgs_ord 1.0000000 0.6812573
## v2cscnsult_ord 0.6812573 1.0000000dataClus=cluster_gabriela[,-1]
row.names(dataClus)=cluster_gabriela$country_namelibrary(cluster)
g.dist = daisy(dataClus, metric="gower")Para ver el número ideal de clusters dependiendo del método, la idea es disminuir dimensiones, si tenemos 10 grupos no se reduciría como si tuviéramos 4 por ejemplo.
En primer lugar para hacer el cálculo de cuantos grupos o cuantos conglomerados se van a constituir recurrí a la medida gap y especifique el método que voy a utilizar y por cada método nos da un número de clusters y a partir de ese número de clusters sugerido por PAM, AGNES Y DIANA, es que los ejecuté y luego se hace la visualización que permite ver cuántos casos están mal agrupado y están quedando fuera de una agrupación correcta. El que mejor agrupa es AGNES.
PAM
# El número sugerido de clusters, específicamente para la estrategia particionante es 9.
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaset.seed(123)
fviz_nbclust(dataClus,
pam,
diss=g.dist,
method = "gap_stat",
k.max = 10,
verbose = F)## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
# AGNES
# El número sugerido de clusters, específicamente para la estrategia AGNES es 9
set.seed(123)
fviz_nbclust(dataClus,
hcut,
diss = g.dist,
method = "gap_stat",
k.max = 10,
verbose = F,
hc_func = "agnes")
# DIANA
# El número sugerido de clusters, específicamente para la estrategia divisiva es 4.
set.seed(123)
fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "diana")
#Evaluar resultados, para saber cuál es la estrategia que más nos
conviene.
###pam
set.seed(123)
res.pam=pam(g.dist,k = 9,cluster.only = F)
###agnes
res.agnes<- hcut(g.dist, k = 9,hc_func='agnes',hc_method = "ward.D")
### diana
res.diana <- hcut(g.dist, k = 4,hc_func='diana')Se escoge el que menos casos abajo porque tienen mal ajuste. Vemos que DIANA tiene más casos abajo
fviz_silhouette(res.pam)## cluster size ave.sil.width
## 1 1 41 0.29
## 2 2 18 0.40
## 3 3 37 0.53
## 4 4 18 0.31
## 5 5 13 0.49
## 6 6 21 0.17
## 7 7 10 0.17
## 8 8 14 0.29
## 9 9 7 0.31
fviz_silhouette(res.agnes)## cluster size ave.sil.width
## 1 1 47 0.33
## 2 2 24 0.32
## 3 3 38 0.53
## 4 4 16 0.33
## 5 5 16 0.35
## 6 6 8 0.29
## 7 7 7 0.61
## 8 8 14 0.17
## 9 9 9 0.36
fviz_silhouette(res.diana)## cluster size ave.sil.width
## 1 1 6 0.37
## 2 2 116 0.26
## 3 3 36 0.45
## 4 4 21 0.30
A partir del gráfico de siluetas, se tiene que el número sugerido para la técnica jerárquica aglomerativa es el que mayor puntuación genera, por lo tanto, es el que menos casos malagrupa.
#Estrategia aglomerativa
set.seed(123)
res.agnes <- hcut(g.dist, k = 9,hc_func='agnes',hc_method = "ward.D")
dataClus$agnes=res.agnes$clusterEs normal que salgan valores negativos porque está estandarizado, se está pidiendo la media de cada grupo por variable.
aggregate(.~ agnes, data=dataClus, mean)## agnes v2elembaut_ord v2mecrit_ord v2dlcountr_ord v2csreprss_ord
## 1 1 0.4662717 -0.06940595 0.1140691 0.4696023
## 2 2 0.5826822 1.10448366 0.4465614 0.5827484
## 3 3 0.9433934 0.89274590 1.0205228 0.9487552
## 4 4 -1.5482422 -2.05645867 -1.5770964 -1.6671174
## 5 5 -1.2026869 -0.90702510 -1.0149693 -1.4114508
## 6 6 -0.2524098 0.24240848 -0.2279912 -0.3376513
## 7 7 -1.5482422 -0.04494992 -0.1155658 -0.2938227
## 8 8 -0.2647511 -0.04494992 -0.3082951 -0.5860131
## 9 9 -1.2410819 -1.19438349 -0.7151681 -1.0989695
## v2cseeorgs_ord v2cscnsult_ord
## 1 0.5180281 -0.03015823
## 2 0.6288955 -0.02883511
## 3 0.8307029 1.52582886
## 4 -1.7798567 -1.36584612
## 5 -1.4357493 -1.45912596
## 6 -0.2559522 -1.45912596
## 7 0.3339463 0.03335145
## 8 -0.6773083 0.03335145
## 9 -1.1517241 0.03335145Se reagrupa y se trata de identificar los grupos mejor posicionados.
dataClus$agnes=dplyr::recode(dataClus$agnes, `3` = 1,
`2`=2,
`1`=3,
`6`=4,
`8`=5,
`5`=6,
`9`=7,
`7`=8,
`4`=9)Proyección de dendograma.
fviz_dend(res.agnes, cex = 0.7, horiz = T)## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
#Estrategia basada en densidad # Como nos interesa clusterizar todos los casos no vamos a usar la estrategia basada en densidad, esta es buena estrategia cuando hay escenarios en los que tenga que excluir los casos y no me conviene dejar fuera 11 cuando con la estrategia AGNES tengo la posibilidad de mal agrupar sólo a 4.
proyeccion = cmdscale(g.dist, k=2 , add = T)dataClus$dim1 <- proyeccion$points[,1]
dataClus$dim2 <- proyeccion$points[,2]base= ggplot(dataClus,aes(x=dim1, y=dim2,label=row.names(dataClus)))
base + geom_text(size=2)
base= ggplot(dataClus,aes(x=dim1, y=dim2)) + coord_fixed()
base + geom_point(size=2, aes(color=as.factor(agnes))) + labs(title = "AGNES")
g.dist.cmd = daisy(dataClus[,c('dim1','dim2')], metric = 'euclidean')library(dbscan)
kNNdistplot(g.dist.cmd, k=7)
library(fpc)##
## Attaching package: 'fpc'## The following object is masked from 'package:dbscan':
##
## dbscandb.cmd = fpc::dbscan(g.dist.cmd, eps=0.065, MinPts=3 , method = 'dist')db.cmd## dbscan Pts=179 MinPts=3 eps=0.065
## 0 1 2 3 4 5 6
## border 11 0 0 0 0 0 0
## seed 0 88 25 27 12 8 8
## total 11 88 25 27 12 8 8dataClus$db=as.factor(db.cmd$cluster)library(ggrepel)
base= ggplot(dataClus[dataClus$db!=0,],aes(x=dim1, y=dim2)) + coord_fixed()
dbplot= base + geom_point(aes(color=db))
dbplot + geom_point(data=dataClus[dataClus$db==0,],
shape=0) 
library(magrittr)#Se observa que son 4 los países mal agrupados.
silPAM=data.frame(res.pam$silinfo$widths)
silPAM$country=row.names(silPAM)
poorPAM=silPAM[silPAM$sil_width<0,'country']%>%sort()
silAGNES=data.frame(res.agnes$silinfo$widths)
silAGNES$country=row.names(silAGNES)
poorAGNES=silAGNES[silAGNES$sil_width<0,'country']%>%sort()
silDIANA=data.frame(res.diana$silinfo$widths)
silDIANA$country=row.names(silDIANA)
poorDIANA=silDIANA[silDIANA$sil_width<0,'country']%>%sort()
library(qpcR)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: minpack.lm## Loading required package: rgl## Loading required package: robustbase## Loading required package: Matrix##
## Attaching package: 'qpcR'## The following object is masked from 'package:DescTools':
##
## RMSEmal_Clus=as.data.frame(qpcR:::cbind.na(poorPAM, poorAGNES,poorDIANA))
mal_Clus## poorPAM poorAGNES poorDIANA
## 1 Australia Bolivia Angola
## 2 Jamaica Lebanon Comoros
## 3 Sri Lanka Mauritania Fiji
## 4 Thailand Tanzania Gabon
## 5 Uzbekistan <NA> Haiti
## 6 Vietnam <NA> India
## 7 <NA> <NA> Jordan
## 8 <NA> <NA> Kuwait
## 9 <NA> <NA> Kyrgyzstan
## 10 <NA> <NA> Maldives
## 11 <NA> <NA> Mali
## 12 <NA> <NA> Mauritania
## 13 <NA> <NA> Morocco
## 14 <NA> <NA> Poland
## 15 <NA> <NA> Serbia# El gráfico muestra de qué manera se agrupo los datos. En algunos casos hay traslape, parece que hay casos encima de otros, esto porque se pide varios grupos para un número no tan alto de casos entonces pareciera que hay un traslape o están ocupando las zonas de otros grupos y como hay 7 variables y 2 dimensiones no necesariamente van a capturar toda la variedad, se reduce a dos solamente para hacer el gráfico.
original=aggregate(.~ agnes, data=dataClus,mean)
proyeccion = cmdscale(g.dist, k=2,add = T)
dataClus$dim1 <- proyeccion$points[,1]
dataClus$dim2 <- proyeccion$points[,2]
base= ggplot(dataClus,aes(x=dim1, y=dim2,label=row.names(dataClus)))
base + geom_text(size=2, aes(color=as.factor(agnes))) + labs(title = "AGNES") 
# II. ANÁLISIS FACTORIAL # Armar la base de apoyo
factor_gabriela = subset(vdem, select = c(country_name, year, v2x_libdem, v2elembaut_ord, v2mecrit_ord, v2dlcountr_ord, v2csreprss_ord, v2cseeorgs_ord, v2cscnsult_ord))factor_gabriela = factor_gabriela[factor_gabriela$year==2021,]factor_gabriela$country_name = NULL
factor_gabriela$year = NULL
factor_gabriela$v2x_libdem = NULLAnálisis Factorial Exploratorio
Explorar las correlaciones entre las variables
# Es otra forma de reducir las dimensiones.
# Las cargas de correlaciones son diferentes a 0.
corMatrix_g = polycor::hetcor(factor_gabriela)$correlations
corMatrix_g## v2elembaut_ord v2mecrit_ord v2dlcountr_ord v2csreprss_ord
## v2elembaut_ord 1.0000000 0.7377899 0.6826106 0.7945379
## v2mecrit_ord 0.7377899 1.0000000 0.6917230 0.7422000
## v2dlcountr_ord 0.6826106 0.6917230 1.0000000 0.7146772
## v2csreprss_ord 0.7945379 0.7422000 0.7146772 1.0000000
## v2cseeorgs_ord 0.7581324 0.7667901 0.6945539 0.8888291
## v2cscnsult_ord 0.6429916 0.6375932 0.6675044 0.7144543
## v2cseeorgs_ord v2cscnsult_ord
## v2elembaut_ord 0.7581324 0.6429916
## v2mecrit_ord 0.7667901 0.6375932
## v2dlcountr_ord 0.6945539 0.6675044
## v2csreprss_ord 0.8888291 0.7144543
## v2cseeorgs_ord 1.0000000 0.6812573
## v2cscnsult_ord 0.6812573 1.0000000Graficar la matriz de correlaciones
# Para saber que las variables covarían se hace una matriz de correlaciones, nos dice que si una aumenta la otra también aumenta o viceversa.Vemos que hay variables que se correlacionan más que otras, pero en general todas se correlacionan bien.
ggcorrplot(corMatrix_g)
Verificar validez del análisis factorial
Verificar si variables se pueden factorizar
Overall MSA es mayor a 0.6, por lo que el análisis factorial es factible.
psych::KMO(corMatrix_g)## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = corMatrix_g)
## Overall MSA = 0.91
## MSA for each item =
## v2elembaut_ord v2mecrit_ord v2dlcountr_ord v2csreprss_ord v2cseeorgs_ord
## 0.93 0.93 0.94 0.85 0.86
## v2cscnsult_ord
## 0.95Descartar una posible matriz de identidad
Sale FALSE (p-value NO es mayor a 0.05), por lo que el análisis factorial es factible.
# Queda comprobado que no es una matriz de identidad, que no se correlacionen entre sí mismas, sino con otras variables además de sí mismas y forman parte de un mismo concepto sin que haya relación de dependencia.
cortest.bartlett(corMatrix_g, n = nrow(factor_gabriela))$p.value>0.05## [1] FALSEDescartar una posible matriz singular
Sale FALSE, por lo que el análisis factorial es factible.
# Queda comprobado que no es una matriz singular.
is.singular.matrix(corMatrix_g)## [1] FALSEDeterminar en cuántos factores se pueden agrupar las variables
# Nos sugiere que el número de factores es 1.
fa.parallel(factor_gabriela, fm = "ML", fa = "fa")
## Parallel analysis suggests that the number of factors = 1 and the number of components = NAObservar las cargas factoriales y ver en qué factores se ubicaría cada variable
#Observamos que el pocertanje de la varianza de las variables es de 72% lo cual está bastante bien. Nos explica la variación de las variables en conjunto, está explicando un 72%.
resfa_g <- fa(factor_gabriela, nfactors = 1, cor = "cor", rotate = "varimax", fm = "minres")
print(resfa_g$loadings, cutoff = 0.5)##
## Loadings:
## MR1
## v2elembaut_ord 0.852
## v2mecrit_ord 0.839
## v2dlcountr_ord 0.801
## v2csreprss_ord 0.925
## v2cseeorgs_ord 0.906
## v2cscnsult_ord 0.772
##
## MR1
## SS loadings 4.347
## Proportion Var 0.724Graficar cómo se agrupan las variables
# Hay un sólo factor que está detrás del comportamiento de esas 6 variables.
fa.diagram(resfa_g)
Evaluar los resultados obtenidos
¿Qué variables aportaron más a los factores?
# Las comunalidades no está diciendo qué porcentaje de las variables con las que estoy trabajando carga bien, la que menor carga es Consulta gubernamental hacia Organizaciones de Sociedad Civil con un 59.65%.
sort(resfa_g$communality)## v2cscnsult_ord v2dlcountr_ord v2mecrit_ord v2elembaut_ord v2cseeorgs_ord
## 0.5965095 0.6420975 0.7045047 0.7260622 0.8212054
## v2csreprss_ord
## 0.8563265Observar los posibles valores proyectados
Para grabar en la base los puntajes de los factores
factor_gabriela$puntaje = resfa_g$scoresAnálisis Factorial Confirmatorio
Construir un modelo lineal
modelog <- ' factorg =~ v2elembaut_ord + v2mecrit_ord + v2dlcountr_ord + v2csreprss_ord + v2cseeorgs_ord + v2cscnsult_ord'Crear un objeto para hacer las validaciones
cfa_fit <- cfa(modelog, data=factor_gabriela,
std.lv=TRUE,
missing="fiml")Preparar los test para las validaciones
allParamCFA=parameterEstimates(cfa_fit,standardized = T)
allFitCFA=as.list(fitMeasures(cfa_fit))Ver si cada variable tiene una buena relación con su factor (p-value menor a 0.05 indica que la variable tiene buena relación con su latente)
allParamCFA[allParamCFA$op=="=~",]## lhs op rhs est se z pvalue ci.lower ci.upper std.lv
## 1 factorg =~ v2elembaut_ord 1.219 0.088 13.906 0 1.048 1.391 1.219
## 2 factorg =~ v2mecrit_ord 0.715 0.054 13.356 0 0.610 0.820 0.715
## 3 factorg =~ v2dlcountr_ord 0.865 0.070 12.307 0 0.727 1.003 0.865
## 4 factorg =~ v2csreprss_ord 1.146 0.069 16.650 0 1.011 1.281 1.146
## 5 factorg =~ v2cseeorgs_ord 1.173 0.072 16.190 0 1.031 1.315 1.173
## 6 factorg =~ v2cscnsult_ord 0.509 0.043 11.907 0 0.425 0.593 0.509
## std.all std.nox
## 1 0.845 0.845
## 2 0.824 0.824
## 3 0.780 0.780
## 4 0.940 0.940
## 5 0.925 0.925
## 6 0.762 0.762Ver si la asignación de variables ha relativamente buena (p-value mayor a 0.05 para validar el modelo)
allFitCFA[c("chisq", "df", "pvalue")]## $chisq
## [1] 29.68009
##
## $df
## [1] 9
##
## $pvalue
## [1] 0.0004972185Otra prueba para ver si el análisis factorial es relativamente bueno (índice Tucker - Lewi debe ser mayor a 0.9)
allFitCFA$tli## [1] 0.9634727Ver si la raíz del error cuadrático medio de aproximación es menor a 0.05 (ver rmsea)
allFitCFA[c('rmsea.ci.lower','rmsea' ,'rmsea.ci.upper')]## $rmsea.ci.lower
## [1] 0.06966505
##
## $rmsea
## [1] 0.1132996
##
## $rmsea.ci.upper
## [1] 0.1597871Hacer predicciones (scores) de las variables latentes
scorescfa=normalize(lavPredict(cfa_fit),
method = "range",
margin=2, # by column
range = c(0, 10))factor_gabriela$prediccion = scorescfaggplot(data=factor_gabriela,aes(x=prediccion,y=puntaje)) + geom_point() + theme_minimal()
cluster_factor_gabriela = merge(dataClus, factor_gabriela$puntaje)