Una de las tareas fundamentales del análisis cuantitativo es garantizar la calidad de las mediciones. Por lo general medimos propiedades o atributos de un constructo. Estos atributos pueden ser de dos tipos:

Variable Latente

En Psicología es muy común que se trabaje con atributos no-observados. Por ejemplo, cuando nos interesamos por la inteligencia, la depresión, la satisfacción, etc.; en todos estos casos tenemos que inferir el valor porque no se observa directamente.

La estrategia más común para medir variables latentes es por medio de escalas: colección de reactivos o ítems cuyas respuestas se combinan en un puntaje que refleja los niveles de la variable latente. El principal objetivo de una escala es discriminar entre individuos según la intensidad de la variable latente.

La teoría clásica de medición parte de una premisa fundamental. A partir del conjunto de reactivos de una escala es posible extraer un puntaje verdadero (true score) que refleja la intensidad de la variable latente. Ahora bien, cada reactivo se compone de dos elementos: puntaje verdadero + error. La clave del análisis factorial es separar ambos componentes.

Veamos el siguiente esquema que resume la teoría de la variable latente:

Noten una consecuencia del diagrama. Lo más común es que nosotros sumamos todas las X para obtener un puntaje total. Ese puntaje total no es el puntaje verdadero. Pero usualmente asumimos que ese puntaje total refleja la variable latente. El problema es que al sumar también estamos sumando el error.

Así es como deberíamos de pensar cada reactivo de una escala; tutorial en spss, por si ocupan. La parte verde es la varianza compartida por todos los reactivos de una escala y, al aislarse, refleja la variable latente. La parte roja es el error, que a veces es varianza compartida por algunos ítems (error sistemático o específico) y no tiene vinculación con la variable latente y, finalmente, está el error aleatorio, que es una particularidad de la persona.

Para poder sumar con tranquilidad una escala necesitamos garantizar unas cuantas cosas, es decir, que no se violen los supuestos. Estas son las condiciones que deben cumplirse:

Muchas veces ignoramos estos problemas. Pero la Psicología se ha destacado por diseñar métodos que nos ayudan a verificar la calidad de nuestras mediciones.

Análisis Factorial Exploratorio

Uno de los primeros análisis que realizamos al analizar una variable latente es el análisis factorial exploratorio, mejor conocido por sus siglas en inglés: EFA.

El Análisis Factorial Exploratorio es un método multivariado que ayuda a examinar las relaciones que ligan a ciertas variables que “deben de ir juntas”. El ejemplo más claro de esto son los ítems de una escala. Cada ítem mide algo diferente, pero juntos miden un fenómeno conjunto: un factor (i.e., la variable latente). El EFA nos ayuda a separar el puntaje verdadero del error y así determinar el número de factores y/o sub-factores que hay en nuestros reactivos. Los principales usos del EFA son:

ENBIARE

Vamos a trabajar con datos de la ENBIARE y analizaremos nuestra vd del semestre: Depresión, con la CESD7. Aquí la variable latente es la depresión. Hasta ahora hemos asumido que funciona bien como escala, pero ahora vamos a explorarlo mediante un análisis factorial exploratorio.

Noten nuestra hipótesis implícita. Estamos diciendo que la CESD7 es unidimensional, lo cual significa que todos los reactivos miden una y sólo una variable latente: depresión. Vamos a corroborarlo.

library(tidyverse)
library(labelled)
library(skimr)
library(janitor)
library(sjPlot)
library(sjmisc)
library(ggrepel)
library(visdat)
library(naniar)
library(psych)
library(GPArotation)
library(lavaan)
library(corrplot)

library(readr)
TSDEM <- read_csv("Datos/TSDEM.csv")

TSDEM1 <- TSDEM %>% 
   mutate(HOGAR= as.character(HOGAR)) %>% 
  mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID


TENBIARE <- read_csv("Datos/TENBIARE.csv")

TENBIARE1 <- TENBIARE %>% 
   mutate(HOGAR= as.character(HOGAR)) %>% 
  mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID

enbiare <- TENBIARE1 %>% # nuevo objeto bd
  left_join(TSDEM1, by= "id_personal")

Para facilitarnos la vida voy a crear un objeto chiquito tipo bd que sólo tenga los reactivos de la cesd7.

cesd <- enbiare %>%
  select(PD2_1 , PD2_2 , PD2_3 , PD2_4 , PD2_5 , PD2_6 , PD2_7)

skim(cesd)
Data summary
Name cesd
Number of rows 31166
Number of columns 7
_______________________
Column type frequency:
numeric 7
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
PD2_1 0 1 0.54 0.84 0 0 0 1 3 ▇▃▁▁▁
PD2_2 0 1 0.61 0.82 0 0 0 1 3 ▇▅▁▂▁
PD2_3 0 1 0.57 0.86 0 0 0 1 3 ▇▃▁▁▁
PD2_4 0 1 0.63 0.91 0 0 0 1 3 ▇▃▁▂▁
PD2_5 0 1 0.92 1.03 0 0 1 2 3 ▇▅▁▃▂
PD2_6 0 1 2.49 0.84 0 2 3 3 3 ▁▁▁▂▇
PD2_7 0 1 0.59 0.85 0 0 0 1 3 ▇▃▁▁▁

Vamos a verlos gráficamente. La idea es que los patrones se vean más o menos igual. Pero lo que buscamos es que las personas que responden 0 en una, respondan 0 o bajos en todas. A la inversa, quien responde 2-3 en una, responda 2-3 en todas. Queremos que haya coherencia en las respuestas.

¿Ven algo raro aquí?

cesd %>% 
  mutate_if(is.numeric, as.factor) %>%
  plot_stackfrq()+
  scale_x_discrete(limits = rev)+
  theme_minimal()

Cuidado al interpretar esta gráfica. Si tenemos un 4.8% de azul profundo en la 1 y 4.0 en la 2, eso no significa que sean las mismas personas quienes respondieron 3 en ambas.

Correlaciones inter-item

Vamos a ver ahora sí estadísticamente qué tenemos. Veamos primero los colores de la siguiente gráfica. Noten el gradiente de color a su mano derecha. Mientras más profundo el azul, más alta la correlación positiva. Mientras más rojo, más alta la correlación negativa.

Al analizar una escala queremos que las respuestas “vayan juntas”, en otras palabras, queremos ver azules profundos porque los ítems suelen codificarse de menor a mayor intensidad en la variable latente.

¿Qué ítem es la que nos preocupa?

cor_pear  <- cor(cesd)
corrplot.mixed(cor_pear, tl.col='black')

Pero aquí el problema también es el tipo de correlación. Vamos a hacer la correlación correcta: policórica. Esta la hacemos cuando tenemos variables de una escala Likert: ordinales. Si tuvieramos solo dicotómicas, usaríamos tetracóricas y si tuvieramos variables con múltiples niveles de medición haríamos una biserial. Para más info de cómo calcularlas, usar el paquete psych.

Noten que, a diferencia de las marices típicas, aquí tenemos unos parámetros adicionales: 3 valores tau. Esta es la frontera en la que, en promedio, los items pasan de 0 a1, de 1 a 2, y de 2 a 3.

cor.policorica  <- cesd %>% 
  polychoric()

cor.policorica
## Call: polychoric(x = .)
## Polychoric correlations 
##       PD2_1 PD2_2 PD2_3 PD2_4 PD2_5 PD2_6 PD2_7
## PD2_1  1.00                                    
## PD2_2  0.74  1.00                              
## PD2_3  0.82  0.74  1.00                        
## PD2_4  0.68  0.70  0.74  1.00                  
## PD2_5  0.44  0.45  0.48  0.44  1.00            
## PD2_6 -0.42 -0.37 -0.44 -0.37 -0.23  1.00      
## PD2_7  0.78  0.68  0.83  0.67  0.50 -0.45  1.00
## 
##  with tau of 
##            1     2     3
## PD2_1  0.366  1.11  1.66
## PD2_2  0.161  1.09  1.75
## PD2_3  0.314  1.07  1.62
## PD2_4  0.257  0.95  1.51
## PD2_5 -0.093  0.62  1.22
## PD2_6 -1.673 -1.11 -0.44
## PD2_7  0.267  1.07  1.64
cor.plot(cor.policorica$rho, numbers = T, upper=FALSE, main = "Correlación inter-item de la CESD-7 en ENBIARE")

Mejoran nuestros números: son más altas las correlaciones cuando usamos la correcta. Pero queda claro que la 6 es problemática. Vamos a seguir para ver cómo se ve este problema con otros indicadores.

KMO

Una de las primeras medidas (en mi opinión demasiado básica) es la medida de Kaiser Meyer Olkin, mejor conocido como KMO. Nos dice si tenemos poder estadístico suficiente para hacer un EFA, en otras palabras, si hay suficiente muestra. En la Enbiare nos sobra. De cualquier manera, el punto de corte es 0.6. Valores menores a 0.6 significaría que debemos quitar ese item o de plano toda la escala. Aquí no es el caso.

KMO(cesd)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cesd)
## Overall MSA =  0.9
## MSA for each item = 
## PD2_1 PD2_2 PD2_3 PD2_4 PD2_5 PD2_6 PD2_7 
##  0.89  0.91  0.86  0.92  0.95  0.96  0.89

ALPHA

Tal vez la medida más conocida en Psicología sea el proverbial alpha. Es una medida mínima para saber la confiabilidad de una prueba. El atractivo es que en una escala de 0-1, nos dice la correlación promedio de los items. El valor aquí es 0.758. Según las interpretaciones clásicas, arriba de 0.7 estamos libres y seguimos y no nos importa nada más. Pero vean cómo hasta el propio R nos da un warning sobre la pregunta 6 (hasta nos dice que debemos recodificarla). Aquí pueden ver qué significa cada estadístico.

alpha_dep <- alpha(cesd)
## Some items ( PD2_6 ) were negatively correlated with the total scale and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
summary(alpha_dep)
## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
##       0.76      0.76    0.81      0.31 3.1 0.0019 0.91 0.56      0.4

Y ahora vean este análisis más específico del alpha. Vean cómo cambia alpha según modificamos un item. En particular miren el raw.r Ese número nos dice la correlación entre ese item y un puntaje sumativo total. Mientras más alto, más relación tiene con depresión. El ítem 3 parece ser especialmente fuerte. Y fuera de la desgracia del 6, el 5 es más o menos débil. También vean r.drop, que sería lo que ganamos o perdemos en alpha al quitar el item.

alpha_dep$item.stats
##           n      raw.r      std.r      r.cor     r.drop      mean        sd
## PD2_1 31166  0.8029827  0.8108612  0.8019992  0.7031984 0.5396265 0.8424133
## PD2_2 31166  0.7773944  0.7858243  0.7517782  0.6713070 0.6132003 0.8193890
## PD2_3 31166  0.8393161  0.8443671  0.8581662  0.7525807 0.5711352 0.8601149
## PD2_4 31166  0.7649547  0.7632443  0.7160820  0.6379030 0.6342809 0.9131644
## PD2_5 31166  0.6257603  0.5979681  0.4652503  0.4239148 0.9169608 1.0284747
## PD2_6 31166 -0.1399583 -0.1324056 -0.3652282 -0.3355769 2.4900212 0.8394722
## PD2_7 31166  0.8022980  0.8060060  0.7977062  0.7003888 0.5877559 0.8544199

OMEGA

Un problema con el alpha es que nos sirve para escalas cuyos items tienen un nivel de medición numérico, pero la mayoría de nuestras escalas tienen un nivel ordinal o dicotómico. Una propuesta es hacer otro tipo de mediciones que se ajustan a la naturaleza de los datos. Con datos ordinales y dicotómicos se recomienda enfocarnos en omega jerárquico. Noten que, de nuevo, mejora nuestras estimaciones y sube a 0.81.

(Favor de ignorar el diagrama)

 omega_dep <- omega(cesd, cor= "poly")
 summary(omega_dep)
## Omega 
## omega(m = cesd, cor = "poly")
## Alpha:                 0.86 
## G.6:                   0.86 
## Omega Hierarchical:    0.81 
## Omega H asymptotic:    0.91 
## Omega Total            0.89 
## 
## With eigenvalues of:
##    g  F1*  F2*  F3* 
## 3.29 0.25 0.36 0.22 
## The degrees of freedom for the model is 3  and the fit was  0 
## The number of observations was  31166  with Chi Square =  61.34  with prob <  0 
## 
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA 
## 
## RMSEA and the  0.9 confidence intervals are  0.025 0.02 0.031
## BIC =  30.3Explained Common Variance of the general factor =  0.8 
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*  F2*  F3*
## Omega total for total scores and subscales    0.89 0.81 0.66 0.89
## Omega general for total scores and subscales  0.81 0.71 0.49 0.69
## Omega group for total scores and subscales    0.06 0.09 0.17 0.21

Lo impresionante es que nos sigue dando valores aceptables, incluso cuando sabemos que hay un ítem malo. Por eso estos índices son criterios laxos.

EFA

El análisis exploratorio se realiza en dos pasos principales – obvio después del análisis exploratorio que acabamos de realizar.

  1. Extracción de factores: determinar el número de variables latentes o subcomponentes que empíricamente mide la escala

  2. Rotación de factores: transformar los resultados para facilitar la interpretación

Extracción de factores

Existen dos métodos de extracción de factores. Unos es el de componentes principales y el otro es el de análisis factorial. El primero se usa para hacer una reducción de dimensiones y sólo obtener un puntaje total. El segundo interesa por medir el tamaño del error y por eso es el de nuestro interés. Esta medición supone separar la varianza en común (factor o variable latente) de la única (error).

Los dos métodos más comunes para extraer los factores son:

  • Máxima verosimilitud - Maximum Likelihood (ML): el más común porque es fácil comparar entre varias combinaciones de ítems con el mismo método de estimación. Es sensible a violaciones de normalidad pero se corrige cuando especificamos la correlación correcta (ej. polychoric).

  • Factor de Ejes principales - Principal Axis Factoring (PAF): preferible con escalas continuas no-normales. Pero menos comparable entre modelos.

Hay varios otros métodos, pero yo les recomiendo quedarse con ML.

Vamos a ver la anatomía del comando:

  • fa= factor analysis

  • 1 = número de factores que vamos a explorar

  • fm= especificar el método

  • ml= Maximum Likelihood

  • cor= tipo de correlación

  • poly= polychoric

fa_cesd <- fa(cesd, 1, fm="ml", cor= "poly")

print(fa_cesd$loadings,cut=.3,digits=2)
## 
## Loadings:
##       ML1  
## PD2_1  0.88
## PD2_2  0.81
## PD2_3  0.93
## PD2_4  0.79
## PD2_5  0.53
## PD2_6 -0.47
## PD2_7  0.88
## 
##                ML1
## SS loadings    4.2
## Proportion Var 0.6

Vamos a interpretar los resultados. Tenemos nuestras 7 preguntas. Cada una ofrece un valor que corresponde al factor “ML1”, esto es, al factor 1 que solicitamos. Para R, el sentido de la variable latente es desconocido. Pero para nosotros, cuando vemos ML1, pensamos que eso significa “Depresión”.

Los valores se llaman cargas factoriales o factor loadings. Este valor va de 0-1 y refleja el cuadrado de la correlación entre el ítem y el puntaje total. En otras palabras, nos dice qué tanta varianza común hay entre el ítem y el factor depresión. Valores más altos indican mayor fuerza. Como ya sabíamos, el ítem 6 no sirve, el 3 es el más fuerte (0.93) y el 5 es el más débil (0.53). Establecí un punto de corte en .3 porque cargas factoriales menores a eso implican que el ítem tiene una relación muy débil con ese factor y no vale la pena considerarlo. Si usan este comando y ven un valor ausente se debe a que su carga es menor a .3.

Implícitamente, al restar 1, cada carga factorial nos da el tamaño del error. Por ejemplo, el error del primer ítem es 1-0.88= 0.12. Por lo tanto, cargas factoriales mayores significan menor error de medición. Por eso nos gustan más esos ítems. si quisiéramos reducir el tamaño de una escala, nos quedaríamos con cargas factoriales altas y quitaríamos las bajas, siempre y cuando se encuentren dentro del mismo factor.

En conjunto, la proporción de varianza compartida de todos los ítems es de 0.6. 60%! Lo que significa que 40% de nuestra varianza es error. Es buena la escala pero en ciencias sociales hay mucha singularidad.

Esta es la forma común de representar el análisis factorial (no rotado); ordenado de mejor ítem a peor.

fa.diagram(fa_cesd, digits=2)

¿Y si no fuera un factor? Vamos a intentar con 2. Noten qué cambia del comando.

fa_cesd2 <- fa(cesd, 2, fm="ml", cor= "poly")

print(fa_cesd2$loadings,cut=.3,digits=2)
## 
## Loadings:
##       ML1   ML2  
## PD2_1  0.79      
## PD2_2        0.98
## PD2_3  0.94      
## PD2_4  0.60      
## PD2_5  0.45      
## PD2_6 -0.51      
## PD2_7  0.96      
## 
##                 ML1  ML2
## SS loadings    3.25 1.04
## Proportion Var 0.46 0.15
## Cumulative Var 0.46 0.61

Aquí estamos “forzando” a que R separe en 2 factores. Pone en un segundo factor el ítem 2, pero solito, lo cual tiene poco sentido; aunque ese ítem explica 15% de la varianza.

Este ejemplo es más o menos fácil porque se trata de un factor, pero no siempre es el caso. Si ustedes anticipan un instrumento con subescalas, tipo Big 5, cada subescala debería acomodarse en un factor distinto (ej. ML1- ML5). Y el objetivo es que su teoría coincida en el acomodo factorial de ítems. Cuando no ocurre, buscan la combinación que funcione mejor para medir diferenciadamente cada subescala.

Veamos un método visual para decidir el número de factores a conservar. No hay un método directo que nos de la respuesta. Así que lo común es calcular varias cosas y ver si convergen. Este comando estima varios métodos y nos da una sugerencia. Aquí recomienda 1 componente, que para nosotros, es equivalente a una variable latente.

fa.parallel(cesd, fm="ml", cor= "poly")

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

Desglosemos un poco esta gráfica. Lo primero es el eje de las X. Ahí nos guía con el número de factores a elegir. En el eje de la Y están esas cosas llamadas eigenvalues. Es un nombre sofisticado que indica la cantidad de varianza que puede explicar un factor. Y nosotros queremos factores que expliquen la mayor cantidad de varianza posible. Esta gráfica indica que el primer factor explica mucha más varianza que el segundo o tercero (cada factor va explicando menos varianza que el precedente). Una regla clásica es que un eigenvalue menor a 1 no vale la pena conservarse porque explica menos que un reactivo solo.

Lo bueno de este comando es que explora varias soluciones para ver si coinciden. Los triángulos reflejan los valores de los reactivos de la enbiare. Los taches usan los valores del fa para simular datos (imaginarios) y ver si coinciden. Las líneas rojas muestran los límites inferiores y superiores de remuestreos tipo bootstrap como los que hacíamos con las pruebas de hipótesis. Afortunadamente, aquí todos coinciden en un factor.

Rotación de factores

La rotación de los factores es una manera de ayudar a la interpretación de los resultados. Es la misma información, pero en términos más comunicables porque pone a todos los factores en los mismos términos. No vamos a entrar a la matemática del asunto, más bien me quiero centrar en una decisión teórica que ayuda a elegir el método de rotación. Pongo en negritas mis recomendaciones de método.

¿Cuando hay dos o más factores, están correlacionados o no? Este supuesto teórico nos lleva a dos familias de métodos:

  • Rotación ortogonal: se asume que los factores son independientes (no están correlacionados). El tipo de rotación sería Varimax, quartimax, equimax, etc.

  • Rotación oblicua: se asume que los factores están correlacionados. El tipo de rotación sería Promax.

Desafortunadamente, mi ejemplo de la cesd-7 no es muy bueno para la rotación porque tiene 1 factor. Pero va igual para que tengan el código. Aunque es algo anti-climático porque es igual que sin rotación : (

fa_cesd_rot <- fa(cesd, 1, fm="ml", cor= "poly", rotate="varimax")

print(fa_cesd_rot$loadings,cut=.3,digits=2)
## 
## Loadings:
##       ML1  
## PD2_1  0.88
## PD2_2  0.81
## PD2_3  0.93
## PD2_4  0.79
## PD2_5  0.53
## PD2_6 -0.47
## PD2_7  0.88
## 
##                ML1
## SS loadings    4.2
## Proportion Var 0.6

Corrección EFA

Ahora sí vamos a hacer bien nuestro análisis factorial y corregimos ese ítem 6 desde el principio. Ahora la bd se llama “dep”.

dep <-  cesd  %>% 
  mutate(PD2_6= case_when(  
    PD2_6== 0 ~ 3,
    PD2_6== 1 ~ 2,
    PD2_6== 2 ~ 1,
    PD2_6== 3 ~ 0))

Frecuencias

Mucho mejor!

dep %>% 
  mutate_if(is.numeric, as.factor) %>%
  plot_stackfrq()+
  scale_x_discrete(limits = rev)+
  theme_minimal()

Correlación iter- item

Vean cómo ahora tenemos todo azul, lo que queremos. Eso no significa que sea maravilloso. Los ítems 5 y 6 se ven pálidos. Eso sugiere que tienen una relación más débil con Depresión.

corrplot.mixed(cor(dep), tl.col='black')

Pero seamos más pro y hagamos las correlaciones poly. Mejor, pero nuestro 6 es un ítem chafita.

cor.policorica2 <- dep %>% 
  polychoric()

cor.plot(cor.policorica2$rho, numbers = T, upper=FALSE, main = "Correlación inter-item de la CESD-7 en ENBIARE")

ALPHA

Sigamos. Ya no hay warning y nos subió a 0.85

alpha_dep1<- alpha(dep)

summary(alpha_dep1)
## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
##       0.85      0.86    0.86      0.46   6 0.0013 0.62 0.64      0.4
alpha_dep1$item.stats
##           n     raw.r     std.r     r.cor    r.drop      mean        sd
## PD2_1 31166 0.8129558 0.8197180 0.8021477 0.7320268 0.5396265 0.8424133
## PD2_2 31166 0.7780528 0.7849122 0.7455420 0.6882218 0.6132003 0.8193890
## PD2_3 31166 0.8531660 0.8576080 0.8591487 0.7854462 0.5711352 0.8601149
## PD2_4 31166 0.7655964 0.7636304 0.7135536 0.6583146 0.6342809 0.9131644
## PD2_5 31166 0.6061072 0.5806428 0.4588346 0.4287716 0.9169608 1.0284747
## PD2_6 31166 0.4959358 0.5045527 0.3609225 0.3355769 0.5099788 0.8394722
## PD2_7 31166 0.8237831 0.8274074 0.8113995 0.7453793 0.5877559 0.8544199

Están más parejos todos. Aunque nos sigue confirmando nuestro diagnóstico de los items fuertes y débiles.

OMEGA

No nos cambia gran cosa; aunque es la forma correcta de medirlo. Recuerden que en otros contextos, con escalas menos probadas o inestables, estos números pueden variar.

 omega_dep1<- omega(dep, cor= "poly")
 summary(omega_dep1)
## Omega 
## omega(m = dep, cor = "poly")
## Alpha:                 0.86 
## G.6:                   0.86 
## Omega Hierarchical:    0.81 
## Omega H asymptotic:    0.91 
## Omega Total            0.89 
## 
## With eigenvalues of:
##    g  F1*  F2*  F3* 
## 3.29 0.25 0.36 0.22 
## The degrees of freedom for the model is 3  and the fit was  0 
## The number of observations was  31166  with Chi Square =  61.34  with prob <  0 
## 
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA 
## 
## RMSEA and the  0.9 confidence intervals are  0.025 0.02 0.031
## BIC =  30.3Explained Common Variance of the general factor =  0.8 
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*  F2*  F3*
## Omega total for total scores and subscales    0.89 0.81 0.66 0.89
## Omega general for total scores and subscales  0.81 0.71 0.49 0.69
## Omega group for total scores and subscales    0.06 0.09 0.17 0.21

EFA

Miren cómo mejoran nuestro números, pero tampoco es dramático. La verdad es que podríamos deshacernos del item 6. ¿Se vale? Esta es una pregunta complicada…

fa_dep <- fa(dep,1, fm="ml", cor= "poly")

print(fa_dep$loadings,cut=.3,digits=2)
## 
## Loadings:
##       ML1 
## PD2_1 0.88
## PD2_2 0.81
## PD2_3 0.93
## PD2_4 0.79
## PD2_5 0.53
## PD2_6 0.47
## PD2_7 0.88
## 
##                ML1
## SS loadings    4.2
## Proportion Var 0.6
fa.diagram(fa_dep, digits=2)

Puntajes de EDA

Una cosa extra que podemos hacer es obtener puntajes. De hecho, una de nuestras motivaciones principales es reducir de 7 números por ítem a 1 resumen que podamos seguir usando con mayor facilidad.

Al inicio hablamos de lo común que es simplemente sumar los ítems y no pensarlo dos veces. Pero ahora ya tenemos más información. Por ejemplo, según el diagrama de arriba, sabemos que el ítem 3 es mucho mejor que el ítem 6 porque nos da más información sobre la depresión de una persona: importa más una respuesta que la otra.

Pero, al sumar, a todas les damos la misma importancia: todas valen igual (0,1,2,3). ¿Cómo usamos mejor esta información? La respuesta es que hacemos un puntaje ponderado. Al igual que con la regresión, usamos las cargas factoriales como parámetros de una ecuación y con eso calculamos un puntaje. Los TAU (fronteras o thresholds) que vimos antes sirven como interceptos y las cargas factoriales como pendientes.

Y, afortunadamente, R lo hace por nosotros.

puntaje <- factor.scores(dep, fa_dep, method = "Thurstone")


puntaje2 <- puntaje$scores

dep2 <- cbind(dep, puntaje2)

summary(dep2$ML1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.9234 -0.8281 -0.3786  0.0000  0.4934  3.7234

Noten que el puntaje factorial está normalizado con media cero.

Vamos a ver cómo se compara con nuestro puntaje original. Rango 0-21 y media 4.3

dep2 <- dep2 %>% 
  mutate(d_total = PD2_1 + PD2_2 + PD2_3 + PD2_4 + PD2_5 + PD2_6 + PD2_7)

summary(dep2$d_total)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   3.000   4.373   7.000  21.000

No es igual, pero sí tiene relación.

dep2 %>% 
  ggplot(aes(x=d_total, y= ML1))+
  geom_point()+
  theme_classic()

Veamos esta gráfica un momento… Esta es la variación del error de medición entre el puntaje sumativo y el puntaje vía el análisis exploratorio.

Seamos más explícitos en nuestras similitudes y diferencias. Vean tanto patrón general como las diferencias particulares. Piensen en particular en la zona del puntaje del 9, que es nuestro punto de corte para determinar si alguien “tiene” o no síntomas de depresión.

dep2 <- dep2 %>%
  mutate(d_total2= ordered(d_total,
                           levels=c("0", "1","2", "3", "4", "5", "6", "7", "8", "9","10","11","12","13","14","15","16","17","18","19","20","21")))


dep2 %>%
  ggplot( aes(x= d_total2, y=ML1)) +
  geom_boxplot()+
    geom_jitter(color="gray", size=0.4, alpha=0.3) +
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

Lo malo con hacer estos puntajes es que… como nos lo advierten los propios creadores de el paquete psych, al crear los scores , este tipo de análisis tiene un problema fundamental. Los resultados te dan una solución estructural de la escala, pero los ponderadores (cargas factoriales) en sí mismos son dependientes de esta muestra. En otras palabras, la enbiare del próximo año dará loadings distintos y puntajes distintos. No es una solución estable como la suma, que al igualar todo, siempre es igual. Por ello el punto de corte se hizo con una suma. Además es más fácil de estimar. El puntaje factorial es más preciso para esta muestra, pero no es generalizable a otras. Además, hay diferentes métodos para estimar estos puntajes, como se explica en el link previo; yo usé “thurstone”, pero hay debate en cuál es el mejor.

De cualquier manera, el EFA nos sirve mucho para conocer nuestros datos y crear y refinar escalas. Y es la antesala para análisis más sofisticados. No se pierdan del análisis confirmatorio (CFA) y de todo el maravilloso mundo de modelar variables latentes (SEM). Para quien quiera entrar, la bibliografía tiene excelentes vías de entrada.

Y, en R, lavaan es el camino.