Evaluación SEM - 2023

Presentación

Estudiante: Corti, Juan Facundo

El presente trabajo constituye la evaluación del curso de posgrado Modelos de ecuaciones estructurales aplicados en la investigación en Psicología. Aplicaciones con R, dictado en la Facultad de Psicología de la Universidad de Buenos Aires. Para su resolución se utilizaron los paquetes readxl, semPlot, MVN y lavaan. Por otro lado, para la construcción de este informe también se utilizaron las librerías knitr, tidyverse, downloadthis, kableExtra y tippy.

Ejercicio 1: Motivos para el consumo de alcohol

Consigna

Un investigador ha desarrollado un cuestionario de 12 ítems (todos ellos tomando valores entre 0 y 8) a los efectos de medir los motivos del consumo de alcohol en adultos jóvenes. El instrumento fue construido para medir tres facetas de este constructo:

  • Motivos de afrontamiento (coping motives): ítems 1 a 4.
  • Motivos sociales (social motives): ítems 5 a 8.
  • Motivos de mejora (enhacement motives): ítems 9 a 12.

Todos los ítems han sido diseñados utilizando frases que puntúan en forma positiva (por ejemplo, el ítem 5: “porque te sientes más seguro de ti mismo”), exceptuando los ítems 11 y 12, que son inversos.

El cuestionario propuesto fue administrado a una muestra de 500 adultos jóvenes.

  1. Realizar el gráfico del modelo propuesto.
  2. Indicar qué parámetros se estimarán en este modelo.
  3. Determinar los grados de libertad del problema.
  4. ¿Es el modelo identificado? Justificar.
  5. Utilizando R, ajustar el modelo propuesto.
  6. Evaluar la significatividad estadística de los parámetros estimados e interpretar.
  7. ¿El estadístico Chi-cuadrado indica un buen ajuste?
  8. Analizar las medidas de bondad de ajuste estudiadas y determinar si se está en presencia de un “buen ajuste”. Justificar la respuesta.
  9. Presentar el gráfico con los estimadores.
  10. Si bien el modelo de tres factores fue diseñado para que éstos no compartan ítems, evaluar si existen relaciones adicionales que mejoren el ajuste.

Base de datos

A continuación se presenta la matriz de varianzas y covarianzas utilizada para la resolución de los ejercicios.

Tabla 1. Matriz de varianzas y covarianzas
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
x1 1.170
x2 0.385 1.071
x3 0.256 0.280 1.073
x4 0.463 0.438 0.462 1.085
x5 0.180 0.252 0.218 0.484 0.934
x6 0.225 0.268 0.268 0.585 0.454 0.940
x7 0.218 0.223 0.251 0.570 0.415 0.481 1.015
x8 0.210 0.250 0.202 0.597 0.470 0.559 0.553 1.105
x9 0.148 0.154 0.114 0.198 0.120 0.130 0.111 0.162 1.045
x10 0.149 0.106 0.065 0.239 0.133 0.188 0.108 0.184 0.472 1.061
x11 0.199 0.139 0.184 0.172 0.044 0.124 0.068 0.066 0.388 0.370 1.054
x12 0.093 0.090 0.104 0.205 0.074 0.171 0.128 0.117 0.362 0.368 0.518 0.991

Resolución

1. Realizar el gráfico del modelo propuesto.

El modelo propuesto es un modelo de medición de la siguiente forma:

mod <- ' afrontamiento =~ x1 + x2 + x3 + x4
          social =~ x5 + x6 + x7 + x8
          mejora =~ x9 + x10 + x11 + x12
        '
fit_ml <- cfa(mod, estimator="ML",
           sample.cov=as.matrix(base), sample.nobs = 500)

semPaths(# Argumentos globales
          fit_ml, what="diagram", layout="tree3", rotation = 4, width=50, height=35, 
         # Etiquetas
          nodeLabels = c(paste0("x",1:12), "Motivos de\nAfrontamiento", "Motivos\nSociales", "Motivos de\nMejora"), label.cex=c(rep(1,12),1,.75,.85), 
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, curvature = 1.5,
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         )

Siendo:

  • \(X_{1, 2, ..., 12}\): las variables observadas.
  • \(\delta_{1, 2, ..., 12}\): sus respectivos errores.
  • Motivos de afrontamiento: la variable latente que predice \(X_{1, 2, 3, 4}\).
  • Motivos sociales: la variable latente que predice \(X_{5, 6, 7, 8}\).
  • Motivos de mejora: la variable latente que predice \(X_{9, 10, 11, 12}\).

2. Indicar qué parámetros se estimarán en este modelo.

Para ajustar este modelo se estimarán los siguientes parámetros:

  • 9 parámetros de regresión (\(\lambda_{1,5,9}\) se fijarán en 1)
  • 12 varianzas de error (una por cada variable observada)
  • 3 varianzas de variables latentes
  • 3 covarianzas entre variables latentes

En total se estimarán 27 parámetros.

3. Determinar los grados de libertad del problema.

Los grados de libertad de un modelo de ecuaciones estructurales están determinados por la diferencia entre la cantidad de momentos diferentes de la matriz de covarianzas de las variables observadas (M; ver fórmula 1) y la cantidad de parámetros a estimar.

\[\begin{equation} \tag{1} M = \frac{p \times (p+1)}{2} \end{equation}\]

Siendo p la cantidad de variables observadas.

De esta manera, el cálculo de los grados de libertad del modelo estará dado por:

\[\begin{align} \tag{2} gl &= \frac{p \times (p+1)}{2} - \text{cant. de parametros} \\ gl &= \frac{12 \times (12+1)}{2} - 27 \\ gl &= \frac{156}{2} - 27 \\ gl &= 78 - 27 \\ gl &= 51 \end{align}\]

4. ¿Es el modelo identificado? Justificar.

Que el modelo esté identificado implica que éste tenga por lo menos tantas ecuaciones como parámetros a estimar (i.e., que tenga grados de libertad positivos). Por ello, el modelo planteado es identificado (\(51\:gl > 0\)).

5. Utilizando R, ajustar el modelo propuesto.

El primer paso para ajustar un modelo es evaluar qué estimador utilizar y, para ello, se debe evaluar la normalidad multivariada. Sin embargo, dado que el dataset utilizado es una matriz de varianzas y covarianzas, y no los datos crudos, no será posible estimarla. Por ello, no se podrían utilizar estimadores que asuman una distribución normal multivariada de los datos (e.g., máxima verosimilitud).

Sin embargo, Gana y Broc ( ) indican que es posible utilizar máxima verosimilitud si las categorías de respuesta superan las 6 opciones. En este caso, las opciones eran 9 (de 0 a 8), por lo que se ajustará el modelo con máxima verosimilitud. De todos modos, también se ajustará el modelo con mínimos cuadrados no ponderados ( ), ya que este estimador produce cargas factoriales más exactas y precisas que mínimos cuadrados ponderados con media y varianza ajustada ( ; ), alternativa frecuentemente sugerida para datos ordinales. Se compararán los indicadores de ajuste para seleccionar aquel que funcione mejor.

mod <- ' afrontamiento =~ x1 + x2 + x3 + x4
          social =~ x5 + x6 + x7 + x8
          mejora =~ x9 + x10 + x11 + x12
        '
fit_ml <- cfa(mod, estimator="ML",
           sample.cov=as.matrix(base), sample.nobs = 500)
fit_uls <- cfa(mod, estimator="ULS",
           sample.cov=as.matrix(base), sample.nobs = 500)



tabla_ajustes <- as.data.frame(rbind(fitmeasures(fit_ml)[c("cfi","tli","rmsea","srmr")],
                    fitmeasures(fit_uls)[c("cfi","tli","rmsea","srmr")]))
tabla_ajustes <- cbind(c("ML","ULS"),tabla_ajustes,
                       c(fitmeasures(fit_ml)[c("chisq")]/fitmeasures(fit_ml)[c("df")],
                         fitmeasures(fit_uls)[c("chisq")]/fitmeasures(fit_uls)[c("df")]))
colnames(tabla_ajustes) <- c("Estimador","CFI","TLI","RMSEA","SRMR","X2/gl")
tabla_ajustes[,2:6] <- round(tabla_ajustes[,2:6],3)

kable(tabla_ajustes,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 2. Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 2. Indicadores de ajuste según estimador
Estimador CFI TLI RMSEA SRMR X2/gl
ML 0.970 0.961 0.044 0.040 1.963
ULS 0.997 0.996 0.019 0.037 1.184

El ajuste con mínimos cuadrados no ponderados fue superior al de máxima verosimilitud en todos los indicadores. En adelante se interpretarán ambos ajustes para evaluar si ambos son válidos.

Ajuste con 

summary(fit_ml, fit.measures = TRUE, standardized=T)
## lavaan 0.6.15 ended normally after 41 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        27
## 
##   Number of observations                           500
## 
## Model Test User Model:
##                                                       
##   Test statistic                               100.117
##   Degrees of freedom                                51
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                              1679.020
##   Degrees of freedom                                66
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.970
##   Tucker-Lewis Index (TLI)                       0.961
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -7845.129
##   Loglikelihood unrestricted model (H1)      -7795.071
##                                                       
##   Akaike (AIC)                               15744.259
##   Bayesian (BIC)                             15858.053
##   Sample-size adjusted Bayesian (SABIC)      15772.354
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.044
##   90 Percent confidence interval - lower         0.031
##   90 Percent confidence interval - upper         0.057
##   P-value H_0: RMSEA <= 0.050                    0.775
##   P-value H_0: RMSEA >= 0.080                    0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.040
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   afrontamiento =~                                                      
##     x1                1.000                               0.472    0.437
##     x2                0.967    0.131    7.369    0.000    0.456    0.441
##     x3                0.994    0.133    7.487    0.000    0.469    0.453
##     x4                2.089    0.223    9.365    0.000    0.985    0.947
##   social =~                                                             
##     x5                1.000                               0.612    0.633
##     x6                1.186    0.091   13.082    0.000    0.725    0.748
##     x7                1.135    0.092   12.361    0.000    0.694    0.690
##     x8                1.251    0.097   12.848    0.000    0.765    0.728
##   mejora =~                                                             
##     x9                1.000                               0.614    0.602
##     x10               1.000    0.107    9.370    0.000    0.614    0.597
##     x11               1.104    0.112    9.874    0.000    0.678    0.661
##     x12               1.076    0.109    9.896    0.000    0.661    0.665
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   afrontamiento ~~                                                      
##     social            0.231    0.033    7.059    0.000    0.802    0.802
##     mejora            0.094    0.020    4.670    0.000    0.325    0.325
##   social ~~                                                             
##     mejora            0.101    0.024    4.167    0.000    0.268    0.268
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1                0.945    0.062   15.332    0.000    0.945    0.809
##    .x2                0.860    0.056   15.315    0.000    0.860    0.805
##    .x3                0.851    0.056   15.271    0.000    0.851    0.795
##    .x4                0.112    0.049    2.264    0.024    0.112    0.103
##    .x5                0.558    0.040   13.805    0.000    0.558    0.599
##    .x6                0.413    0.035   11.941    0.000    0.413    0.440
##    .x7                0.531    0.041   13.071    0.000    0.531    0.524
##    .x8                0.518    0.042   12.382    0.000    0.518    0.469
##    .x9                0.665    0.053   12.483    0.000    0.665    0.638
##    .x10               0.681    0.054   12.567    0.000    0.681    0.644
##    .x11               0.592    0.053   11.215    0.000    0.592    0.563
##    .x12               0.551    0.050   11.121    0.000    0.551    0.558
##     afrontamiento     0.223    0.046    4.848    0.000    1.000    1.000
##     social            0.374    0.051    7.308    0.000    1.000    1.000
##     mejora            0.378    0.060    6.250    0.000    1.000    1.000

Ajuste con 

summary(fit_uls, fit.measures = TRUE, standardized=T)
## lavaan 0.6.15 ended normally after 41 iterations
## 
##   Estimator                                        ULS
##   Optimization method                           NLMINB
##   Number of model parameters                        27
## 
##   Number of observations                           500
## 
## Model Test User Model:
##                                                       
##   Test statistic                                60.379
##   Degrees of freedom                                51
##   P-value (Unknown)                                 NA
## 
## Model Test Baseline Model:
## 
##   Test statistic                              2971.009
##   Degrees of freedom                                66
##   P-value                                           NA
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997
##   Tucker-Lewis Index (TLI)                       0.996
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.019
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.036
##   P-value H_0: RMSEA <= 0.050                    1.000
##   P-value H_0: RMSEA >= 0.080                    0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.037
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   afrontamiento =~                                                      
##     x1                1.000                               0.485    0.448
##     x2                1.034    0.083   12.400    0.000    0.501    0.485
##     x3                0.978    0.081   12.120    0.000    0.474    0.457
##     x4                2.071    0.149   13.892    0.000    1.004    0.964
##   social =~                                                             
##     x5                1.000                               0.614    0.636
##     x6                1.208    0.078   15.461    0.000    0.742    0.765
##     x7                1.126    0.074   15.216    0.000    0.692    0.687
##     x8                1.223    0.079   15.500    0.000    0.752    0.715
##   mejora =~                                                             
##     x9                1.000                               0.643    0.629
##     x10               1.010    0.091   11.075    0.000    0.649    0.630
##     x11               1.000    0.090   11.060    0.000    0.643    0.626
##     x12               0.988    0.090   11.038    0.000    0.635    0.638
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   afrontamiento ~~                                                      
##     social            0.222    0.018   12.657    0.000    0.746    0.746
##     mejora            0.111    0.012    9.167    0.000    0.356    0.356
##   social ~~                                                             
##     mejora            0.107    0.012    8.629    0.000    0.270    0.270
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1                0.935    0.052   17.894    0.000    0.935    0.799
##    .x2                0.819    0.053   15.512    0.000    0.819    0.765
##    .x3                0.849    0.052   16.354    0.000    0.849    0.791
##    .x4                0.077    0.092    0.836    0.403    0.077    0.071
##    .x5                0.557    0.057    9.756    0.000    0.557    0.596
##    .x6                0.390    0.064    6.112    0.000    0.390    0.415
##    .x7                0.536    0.061    8.806    0.000    0.536    0.528
##    .x8                0.540    0.064    8.390    0.000    0.540    0.489
##    .x9                0.632    0.067    9.489    0.000    0.632    0.605
##    .x10               0.639    0.067    9.529    0.000    0.639    0.603
##    .x11               0.640    0.067    9.613    0.000    0.640    0.608
##    .x12               0.587    0.066    8.906    0.000    0.587    0.593
##     afrontamiento     0.235    0.027    8.725    0.000    1.000    1.000
##     social            0.377    0.035   10.669    0.000    1.000    1.000
##     mejora            0.413    0.049    8.386    0.000    1.000    1.000

6. Evaluar la significatividad estadística de los parámetros estimados e interpretar.

Todos los parámetros de regresión fueron significativos y, en su versión estandarizada, se encuentran dentro de los . Sin embargo, en la consigna se establece que los ítems 11 y 12 son inversos, por lo que las estimaciones de los parámetros estandarizados deberían estar entre -1 y 0, pero ambos son positivos (0.661 y 0.665, respectivamente).

Las varianzas y covarianzas estimadas también resultaron todas significativas y positivas.

Los parámetros de regresión también fueron significativos con este estimador y dentro de . Los ítems 11 y 12 también obtuvieron estimaciones estandarizadas positivas (0.626 y 0.638, respectivamente).

Las varianzas y covarianzas fueron todas significativamente distintas de 0, a excepción de la varianza de x4 (p = 0.403). Esto es consistente con la estimación del parámetro de regresión estandarizado para esa variable (0.964), que implicaría una correlación casi perfecta con el constructo Motivos de Afrontamiento.

7. ¿El estadístico Chi-cuadrado indica un buen ajuste?

El p-valor asociado a un estadístico \(\chi^2\) de 100.117 con 51 grados de libertad es menor a 0.001, por lo que no indicaría un buen ajuste.

El estadístico obtenido para el test de \(\chi^2\) con fue 60.379. Teniendo 51 grados de libertad, la probabilidad de obtener ese nivel o uno mayor es . Que este test no sea significativo es indicador de un buen ajuste del modelo.

8. Analizar las medidas de bondad de ajuste estudiadas y determinar si se está en presencia de un “buen ajuste”. Justificar la respuesta.

En la Tabla 2 se presentaron los principales indicadores de bondad de ajuste para . De acuerdo con Gana y Broc ( ), ambos ajustes presentan muy buenos indicadores de bondad de ajuste absolutos (\(SRMR \le 0.08\)), parsimoniosos (\(RMSEA \le 0.05\)) e incrementales (\(CFI \ge 0.95\) y \(TLI \ge 0.95\)).

9. Presentar el gráfico con los estimadores.

A continuación se presenta el modelo con las estimaciones estandarizadas obtenidas tanto con como con .

semPaths(# Argumentos globales
          fit_ml, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          nodeLabels = c(paste0("x",1:12), "Motivos de\nAfrontamiento", "Motivos\nSociales", "Motivos de\nMejora"), label.cex=c(rep(1,12),1,.75,.85), edge.label.cex = c(rep(1,12),rep(.75,12),rep(1,6)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, curvature = 1.5,
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

plot <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          nodeLabels = c(paste0("x",1:12), "Motivos de\nAfrontamiento", "Motivos\nSociales", "Motivos de\nMejora"), label.cex=c(rep(1,12),1,.75,.85), edge.label.cex = c(rep(1,12),rep(.75,12),rep(1,6)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, curvature = 1.5,
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         )

10. Si bien el modelo de tres factores fue diseñado para que éstos no compartan ítems, evaluar si existen relaciones adicionales que mejoren el ajuste.

Para evaluar si el agregado de nuevas relaciones mejoran el ajuste se pueden indagar los índices de modificación, considerando que valores superiores a 5 en estos índices pueden ser de relevancia para el modelo.

##        lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
## 120    x11 ~~ x12 26.755  0.255   0.255    0.446    0.446
## 42  social =~  x4 24.509  1.935   1.183    1.137    1.137
## 115     x9 ~~ x10 19.886  0.194   0.194    0.287    0.287
## 55      x1 ~~  x2 19.091  0.184   0.184    0.204    0.204
## 67      x2 ~~  x4  8.555 -0.128  -0.128   -0.414   -0.414
## 39  social =~  x1  8.296 -0.543  -0.332   -0.307   -0.307
## 83      x3 ~~ x11  8.153  0.103   0.103    0.145    0.145
## 117     x9 ~~ x12  7.684 -0.125  -0.125   -0.206   -0.206
## 118    x10 ~~ x11  7.504 -0.127  -0.127   -0.200   -0.200
## 64      x1 ~~ x11  6.543  0.097   0.097    0.130    0.130
## 80      x3 ~~  x8  5.404 -0.078  -0.078   -0.118   -0.118
## 119    x10 ~~ x12  5.371 -0.105  -0.105   -0.171   -0.171

##        lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
## 42  social =~  x4 18.328  1.369   0.841    0.807    0.807
## 55      x1 ~~  x2 12.281  0.174   0.174    0.198    0.198
## 120    x11 ~~ x12 12.157  0.222   0.222    0.363    0.363
## 39  social =~  x1  8.956 -0.525  -0.323   -0.298   -0.298

En ambos ajustes los índices de modificación sugieren agregar términos de regresión entre Motivos sociales y X1 y X4, y covariancias entre X1 y X2 y entre X11 y X12. En el ajuste obtenido con máxima verosimilitud, además de esos términos, también se identificaron otras 8 covarianzas que podrían mejorar el ajuste.

Más allá del valor del índice de modificación del ajuste, es relevante evaluar si los términos incluidos son significativos.

Términos de regresión

Primero se probaron los términos de regresión entre Motivos sociales y X1 y X4.

mod2 <- ' afrontamiento =~ x1 + x2 + x3 + x4
          social =~ x5 + x6 + x7 + x8 + x1 + x4
          mejora =~ x9 + x10 + x11 + x12
        '
fit_ml2 <- cfa(mod2, estimator="ML",
           sample.cov=as.matrix(base), sample.nobs = 500)
fit_uls2 <- cfa(mod2, estimator="ULS",
           sample.cov=as.matrix(base), sample.nobs = 500)

tab_regr_ml <- summary(fit_ml2, standardized=T)$pe[summary(fit_ml2)$pe$lhs=="social"&
                                    (summary(fit_ml2)$pe$rhs=="x1"|
                                       summary(fit_ml2)$pe$rhs=="x4"),c(3,5:10)]
tab_regr_ml[2:7] <- round(tab_regr_ml[2:7],3)
tab_regr_ml <- cbind("modelo"=rep("ML",2),tab_regr_ml)

tab_regr_uls <- summary(fit_uls2, standardized=T)$pe[summary(fit_uls2)$pe$lhs=="social"&
                                    (summary(fit_uls2)$pe$rhs=="x1"|
                                       summary(fit_uls2)$pe$rhs=="x4"),c(3,5:10)]
tab_regr_uls[2:7] <- round(tab_regr_uls[2:7],3)
tab_regr_uls <- cbind("modelo"=rep("ULS",2),tab_regr_uls)

tabla_regr <- rbind(tab_regr_ml,tab_regr_uls)
colnames(tabla_regr)[2] <- "termino"
tabla_regr$termino <- rep(c("social =~ x1", "social =~ x4"),2)

kable(tabla_regr, row.names = F,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 3. Términos de regresión agregados según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 3. Términos de regresión agregados según estimador
modelo termino est se z pvalue std.lv std.all
ML social =~ x1 -0.265 0.178 -1.488 0.137 -0.162 -0.150
ML social =~ x4 0.723 0.125 5.786 0.000 0.442 0.424
ULS social =~ x1 -0.304 0.213 -1.423 0.155 -0.186 -0.172
ULS social =~ x4 0.711 0.154 4.613 0.000 0.437 0.419

En ambos escenarios, sólo la relación entre Motivos sociales y x4 resultó significativa.

Términos de covarianza

La inclusión de estos términos se realizó sin considerar la relación entre Motivos sociales y x4, sino tomando de base el modelo original.

mod3 <- ' afrontamiento =~ x1 + x2 + x3 + x4
          social =~ x5 + x6 + x7 + x8
          mejora =~ x9 + x10 + x11 + x12
          x1 ~~  x2
          x1 ~~ x11
          x2 ~~  x4
          x3 ~~  x8
          x3 ~~ x11
          x9 ~~ x10
          x9 ~~ x12
          x10 ~~ x11
          x10 ~~ x12
          x11 ~~ x12
        '
mod4 <- ' afrontamiento =~ x1 + x2 + x3 + x4
          social =~ x5 + x6 + x7 + x8
          mejora =~ x9 + x10 + x11 + x12
          x1 ~~  x2
          x11 ~~ x12
        '
fit_ml3 <- cfa(mod3, estimator="ML",
           sample.cov=as.matrix(base), sample.nobs = 500)
fit_uls3 <- cfa(mod4, estimator="ULS",
           sample.cov=as.matrix(base), sample.nobs = 500)

tab_cov_ml <- summary(fit_ml3, standardized=T)$pe[13:22,c(1:3,5:10)]
tab_cov_ml[4:9] <- round(tab_cov_ml[4:9],3)
tab_cov_ml$termino <- paste0(tab_cov_ml$lhs,tab_cov_ml$op,tab_cov_ml$rhs)
tab_cov_ml$modelo <- rep("ML", nrow(tab_cov_ml))
tab_cov_ml <- tab_cov_ml[c(11,10,4:9)]

tab_cov_uls <- summary(fit_uls3, standardized=T)$pe[13:14,c(1:3,5:10)]
tab_cov_uls[4:9] <- round(tab_cov_uls[4:9],3)
tab_cov_uls$termino <- paste0(tab_cov_uls$lhs,tab_cov_uls$op,tab_cov_uls$rhs)
tab_cov_uls$modelo <- rep("ULS", nrow(tab_cov_uls))
tab_cov_uls <- tab_cov_uls[c(11,10,4:9)]

tabla_cov <- rbind(tab_cov_ml,tab_cov_uls)

kable(tabla_cov, row.names = F,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 4. Términos de covarianza agregados según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 4. Términos de covarianza agregados según estimador
modelo termino est se z pvalue std.lv std.all
ML x1~~x2 0.158 0.045 3.491 0.000 0.158 0.177
ML x1~~x11 0.091 0.037 2.492 0.013 0.091 0.108
ML x2~~x4 -0.069 0.047 -1.471 0.141 -0.069 -0.471
ML x3~~x8 -0.059 0.034 -1.755 0.079 -0.059 -0.089
ML x3~~x11 0.094 0.035 2.663 0.008 0.094 0.118
ML x9~~x10 -0.078 0.155 -0.504 0.614 -0.078 -0.160
ML x9~~x12 -0.105 0.120 -0.875 0.382 -0.105 -0.191
ML x10~~x11 -0.078 0.106 -0.739 0.460 -0.078 -0.141
ML x10~~x12 -0.180 0.186 -0.967 0.334 -0.180 -0.386
ML x11~~x12 0.136 0.102 1.331 0.183 0.136 0.219
ULS x1~~x2 0.173 0.049 3.541 0.000 0.173 0.191
ULS x11~~x12 0.211 0.057 3.675 0.000 0.211 0.295

En este caso, los resultados varían un poco más entre ambos ajustes. Si bien la covarianza entre x1 y x2 fue significativa en ambos modelos, la covarianza entre x11 y x12 sólo fue significativa en el modelo ajustado por mínimos cuadrados no ponderados. En el modelo de máxima verosimilitud también fueron significativas las covariancias entre x1 y x11 y x3 y x11.

Conclusiones

En términos de las medidas de bondad de ajuste, tanto el modelo ajustado con como el ajustado con presentan valores óptimos, beneficiando levemente al ajustado con . Contrariamente, el análisis de los coeficientes estimados muestra que este último presenta una varianza no significativa en x4, a diferencia del ajustado con máxima verosimilitud, cuyos coeficientes son todos significativos. Por otro lado, ambos ajustes estimaron coeficientes de regresión positivos para x11 y x12, resultado que contradice el modelo teórico, dado que en el instrumento son ítems inversos. Sin embargo, esto puede deberse a una corrección realizada en la codificación al cargar la base original y no a una falla del instrumento.

Considerando los índices de modificación del ajuste, podrían incluirse algunos términos al modelo para mejorarlo. Sin embargo, teniendo en cuenta que esto entraría en conflicto con el modelo teórico del que se parte y, además, que el modelo original tiene muy buenos indicadores de ajuste, no se considera recomendable este paso.

Por todo esto, se recomienda considerar como modelo final al modelo original ajustado por máxima verosimilitud. No obstante, cabe recordar que no fue corroborado el supuesto de distribución normal multivariada que asume este estimador, por lo que se sugiere su análisis con la base de datos cruda antes de reportar resultados.

Ejercicio 2: Democracia política e industrialización

Consigna

Se cuenta con un archivo de datos que contiene variables de democracia política e industrialización en 75 países en desarrollo de la década de 1960. Los datos están disponibles en el paquete Lavaan (Data “PoliticalDemocracy”)

Las variables contenidas en dicho dataframe son las siguientes:

  • y1: Valoración de expertos sobre la libertad de prensa en 1960.
  • y2: La libertad de la oposición política en 1960.
  • y3: La equidad de las elecciones en 1960.
  • y4: La eficacia de la legislatura elegida en 1960.
  • y5: Valoración de expertos sobre la libertad de prensa en 1965.
  • y6: La libertad de la oposición política en 1965.
  • y7: La equidad de las elecciones en 1965
  • y8: La eficacia de la legislatura elegida en 1965.
  • x1: PBI per capita in 1960.
  • x2: Consumo de energía inanimada per cápita en 1960.
  • x3: El porcentaje de la población económicamente activa en la industria en 1960.

Los investigadores proponen el siguiente modelo:

  1. Analizar e interpretar las relaciones propuestas por los investigadores
  2. Escribir las ecuaciones del modelo.
  3. Utilizando la función sem() correr este modelo de ecuaciones estructurales usando Máxima Verosimilitud.
  4. Analizar los parametros obtenidos
  5. Utilizando la función semPaths() realizar el gráfico del modelo en R.
  6. Evaluar la bondad de ajuste del modelo
  7. Evaluar la normalidad de las variables observables.
  8. Utilizar un método de estimación adecuado en caso de no cumplirse el supuesto de normalidad multivariada.

(Sug. Utilizar el script “curso_SEM_c3.R” como guía general para la resolución del problema)

Base de datos

A continuación se presenta la base de datos utilizada para la resolución de los ejercicios.

Resolución

1. Analizar e interpretar las relaciones propuestas por los investigadores

El modelo propuesto se puede descomponer en dos partes: los modelos de medición y el modelo estructural.

Los modelos de medición están compuestos por aquellas ecuaciones que permiten inferir un factor latente (no observable) a partir de variables manifiestas (observables). En este caso, se presentan tres modelos de medición: uno que supone que los valores de \(y_1\) a \(y_4\) se pueden predecir a partir de un factor latente denominado \(dem60\), otro que propone a la variable latente \(dem65\) como predictora de los valores de \(y_5\) a \(y_8\), y finalmente un factor inobservable (\(ind60\)) para predecir las variables manifiestas \(x_1\) a \(x_3\).

Por su parte, los modelos estructurales establecen relaciones entre las variables latentes. Los investigadores de este estudio, por ejemplo, establecen que los niveles de \(ind60\) permiten predecir los niveles de \(dem60\) y que estas dos variables, en conjunto, predicirían los niveles de \(dem65\).

Para ver la representación matemática de estas escuaciones ver el punto 2 de este ejercicio.

Es importante considerar que las variables \(y_5\) a \(y_8\) son mediciones 5 años después de las variables \(y_1\) a \(y_4\). Esto conlleva una falta de independencia entre esas observaciones que es necesario modelar al establecer las relaciones a incluir en el modelo. Esto se logra al incorporar en el modelo la covariación entre estas variables, o correlaciones entre sus residuos (flechas que relacionan \(y_1\) con \(y_5\), \(y_2\) con \(y_6\), etc.).

De esta manera, la sintaxis para ajustar el modelo con lavaan sería:

modelo <- '
            # Ecuaciones de medicion
             dem60 =~ y1 + y2 + y3 + y4
             dem65 =~ y5 + y6 + y7 + y8
             ind60 =~ x1 + x2 + x3
            
            # Ecuaciones estructurales
             dem60 ~ ind60
             dem65 ~ ind60 + dem60 
             
            # Correlaciones de residuos
             y1 ~~ y5
             y2 ~~ y6
             y3 ~~ y7
             y4 ~~ y8
                                          '

2. Escribir las ecuaciones del modelo.

En los modelos , las ecuaciones que se proponen para el modelo a evaluar pueden dividirse en ecuaciones de medición (aquellas que vinculan las variables observables con las variables latentes) y ecuaciones estructurales (aquellas que vinculan variables latentes entre sí).

Ecuaciones de medición

Dem60

\[\begin{equation} y_1 = \lambda_1 \eta_1 + \delta_1 \\ y_2 = \lambda_2 \eta_1 + \delta_2 \\ y_3 = \lambda_3 \eta_1 + \delta_3 \\ y_4 = \lambda_4 \eta_1 + \delta_4 \end{equation}\]

Siendo:

  • \(y_1\) a \(y_4\) las variables observables.
  • \(\lambda_1\) a \(\lambda_4\) los coeficientes de regresión.
  • \(\eta_1\) la variable latente endógena Dem60.
  • \(\delta_1\) a \(\delta_4\) los errores.
Dem65

\[\begin{equation} y_5 = \lambda_5 \eta_2 + \delta_5 \\ y_6 = \lambda_6 \eta_2 + \delta_6 \\ y_7 = \lambda_7 \eta_2 + \delta_7 \\ y_8 = \lambda_8 \eta_2 + \delta_8 \end{equation}\]

Siendo:

  • \(y_5\) a \(y_8\) las variables observables.
  • \(\lambda_5\) a \(\lambda_8\) los coeficientes de regresión.
  • \(\eta_2\) la variable latente endógena Dem65.
  • \(\delta_5\) a \(\delta_8\) los errores.
Ind60

\[\begin{equation} x_1 = \lambda_9 \xi_1 + \delta_9 \\ x_2 = \lambda_{10} \xi_1 + \delta_{10} \\ x_3 = \lambda_{11} \xi_1 + \delta_{11} \end{equation}\]

Siendo:

  • \(x_1\) a \(x_3\) las variables observables.
  • \(\lambda_9\) a \(\lambda_{11}\) los coeficientes de regresión.
  • \(\xi_1\) la variable latente exógena Ind60.
  • \(\delta_9\) a \(\delta_{11}\) los errores.

Ecuaciones estructurales

\[\begin{equation} \eta_1 = \gamma_1 \xi_1 + \zeta_1 \\ \eta_2 = \gamma_2 \xi_1 + \zeta_2 \end{equation}\]

Siendo:

  • \(\eta_1\) y \(\eta_2\) las variables latentes endógenas.
  • \(\gamma_1\) a \(\gamma_2\) los coeficientes de regresión.
  • \(\xi_1\) la variable latente exógena Ind60.
  • \(\zeta_1\) a \(\zeta_2\) los errores.

3. Utilizando la función sem() correr este modelo de ecuaciones estructurales usando Máxima Verosimilitud.

mod <- '# Ecuaciones de medicion
             dem60 =~ y1 + y2 + y3 + y4
             dem65 =~ y5 + y6 + y7 + y8
             ind60 =~ x1 + x2 + x3
            
        # Ecuaciones estructurales
             dem60 ~ ind60
             dem65 ~ ind60 + dem60 
             
        # Correlaciones de residuos
             y1 ~~ y5
             y2 ~~ y6
             y3 ~~ y7
             y4 ~~ y8
'

# Ajuste del modelo con ML
fit <- sem(mod, data = base2)
summary(fit, fit.measures = T, standardized=T)
## lavaan 0.6.15 ended normally after 58 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        29
## 
##   Number of observations                            75
## 
## Model Test User Model:
##                                                       
##   Test statistic                                50.835
##   Degrees of freedom                                37
##   P-value (Chi-square)                           0.064
## 
## Model Test Baseline Model:
## 
##   Test statistic                               730.654
##   Degrees of freedom                                55
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.980
##   Tucker-Lewis Index (TLI)                       0.970
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1554.146
##   Loglikelihood unrestricted model (H1)      -1528.728
##                                                       
##   Akaike (AIC)                                3166.292
##   Bayesian (BIC)                              3233.499
##   Sample-size adjusted Bayesian (SABIC)       3142.099
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.071
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.115
##   P-value H_0: RMSEA <= 0.050                    0.234
##   P-value H_0: RMSEA >= 0.080                    0.396
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.050
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   dem60 =~                                                              
##     y1                1.000                               2.145    0.824
##     y2                1.388    0.188    7.401    0.000    2.977    0.760
##     y3                1.053    0.161    6.552    0.000    2.259    0.694
##     y4                1.368    0.153    8.928    0.000    2.933    0.881
##   dem65 =~                                                              
##     y5                1.000                               2.014    0.777
##     y6                1.317    0.180    7.314    0.000    2.654    0.790
##     y7                1.326    0.174    7.618    0.000    2.672    0.817
##     y8                1.391    0.171    8.118    0.000    2.803    0.870
##   ind60 =~                                                              
##     x1                1.000                               0.670    0.920
##     x2                2.181    0.139   15.720    0.000    1.460    0.973
##     x3                1.819    0.152   11.966    0.000    1.218    0.872
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   dem60 ~                                                               
##     ind60             1.435    0.385    3.728    0.000    0.448    0.448
##   dem65 ~                                                               
##     ind60             0.507    0.209    2.425    0.015    0.168    0.168
##     dem60             0.816    0.100    8.156    0.000    0.869    0.869
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .y1 ~~                                                                 
##    .y5                0.892    0.366    2.433    0.015    0.892    0.370
##  .y2 ~~                                                                 
##    .y6                1.893    0.762    2.486    0.013    1.893    0.361
##  .y3 ~~                                                                 
##    .y7                1.268    0.623    2.035    0.042    1.268    0.287
##  .y4 ~~                                                                 
##    .y8                0.141    0.464    0.303    0.762    0.141    0.056
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .y1                2.181    0.456    4.779    0.000    2.181    0.322
##    .y2                6.490    1.231    5.271    0.000    6.490    0.423
##    .y3                5.490    0.991    5.538    0.000    5.490    0.518
##    .y4                2.470    0.660    3.741    0.000    2.470    0.223
##    .y5                2.662    0.506    5.260    0.000    2.662    0.396
##    .y6                4.249    0.817    5.201    0.000    4.249    0.376
##    .y7                3.560    0.712    4.999    0.000    3.560    0.333
##    .y8                2.531    0.609    4.159    0.000    2.531    0.244
##    .x1                0.082    0.020    4.177    0.000    0.082    0.154
##    .x2                0.120    0.070    1.708    0.088    0.120    0.053
##    .x3                0.467    0.090    5.172    0.000    0.467    0.239
##    .dem60             3.678    0.885    4.154    0.000    0.799    0.799
##    .dem65             0.350    0.187    1.865    0.062    0.086    0.086
##     ind60             0.448    0.087    5.171    0.000    1.000    1.000

4. Analizar los parametros obtenidos

Todas las estimaciones de los parámetros de regresión obtenidas fueron significativas y sus versiones estandarizadas se encuentran dentro del rango aceptable (0 - 1). En los modelos de medición, las correlaciones de los ítems con los factores latentes es alta. Se destaca el coeficiente estandarizado de \(x_2\), con una correlación casi perfecta, lo que podría ser un problema. En el modelo estructural hay mayor variedad en la intensidad de las relaciones, pero aun así son todas significativas.

En lo que respecta a las varianzas estimadas de las variables observables, sólo la de \(x_2\) no fue significativa. Considerando la fuerte correlación entre esta variable y el factor latente \(ind60\), era de esperar que su varianza no difiera significativamente de cero. En cuanto a los factores, la varianza de \(dem65\) tampoco resultó significativa.

Finalmente, del análisis de las covarianzas de los ítems se desprende que la correlación de los residuos de \(x_4\) y \(x_8\) no fue significativa. El resto de las covarianzas sí lo fueron.

5. Utilizando la función semPaths() realizar el gráfico del modelo en R.

semPaths(# Argumentos globales
          fit, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,13),.85), 
          edge.label.cex = c(rep(.92,11),rep(1.2,3),rep(.87,4),rep(.75,17)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, curvature = c(1.85, 2.45, 3.05, 3.65, 1.25),
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

6. Evaluar la bondad de ajuste del modelo

En el punto 3 se presentaron los principales indicadores de bondad de ajuste para el modelo estimado con máxima verosimilitud. De acuerdo con Gana y Broc ( ), el ajuste presenta muy buenos indicadores de bondad de ajuste absolutos (\(SRMR \le 0.08\)) e incrementales (\(CFI \ge 0.95\) y \(TLI \ge 0.95\)). En cuanto a indicadores de bondad de ajuste parsimoniosos (\(RMSEA \le 0.05\)), el ajuste es bueno aunque marginal (\(RMSEA = 0.071\)).

7. Evaluar la normalidad de las variables observables.

normalidad <- mvn(data = base2,mvnTest = "mardia",
              univariateTest = "SW",univariatePlot = "qqplot",
              multivariatePlot = "qq")

En el gráfico de las distancias de Mahalanobis al cuadrado vs. los cuantiles teóricos de una distribución \(\chi^2\) se pueden ver indicios de incumplimiento del supuesto de normalidad multivariada que se asume al utilizar máxima verosimilitud como estimador. Se aprecia cierta estructura en la distribución de puntos que se aleja de la recta esperada si los datos se distribuyeran de forma normal multivariada, probablemente debido a la naturaleza ordinal de los datos originales.

Por otro lado, en los gráficos de cuantines teóricos de una distribución normal vs. los cuantiles empíricos se observan grandes desviaciones de las correspondientes distribuciones normales univariadas en las variables \(y_1\) a \(y_8\). Los gráficos de las variables \(x_1\) a \(x_3\) no parecieran presentar problemas graves en torno a la normalidad de su distribución.

tabla <- as.data.frame(normalidad$multivariateNormality)
tabla$Statistic <- c(round(as.numeric(levels(tabla$Statistic)),2),NA)
tabla$`p value` <- c(round(as.numeric(levels(tabla$`p value`)),3),NA)

kable(tabla,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 5. Indicadores de normalidad multivariada") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 5. Indicadores de normalidad multivariada
Test Statistic p value Result
Mardia Skewness 344.49 0.010 NO
Mardia Kurtosis -1.22 0.222 YES
MVN NO

Coincidiendo con los gráficos, el test de hipótesis de Mardia para evaluar la asimetría resultó significativo, por lo que hay evidencias de incumplimiento del supuesto de normalidad multivariada.

kable(normalidad$univariateNormality,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 6. Indicadores de normalidad univariada") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 6. Indicadores de normalidad univariada
Test Variable Statistic p value Normality
Shapiro-Wilk y1 0.9316 6e-04 NO
Shapiro-Wilk y2 0.8311 <0.001 NO
Shapiro-Wilk y3 0.8547 <0.001 NO
Shapiro-Wilk y4 0.8994 <0.001 NO
Shapiro-Wilk y5 0.9598 0.0179 NO
Shapiro-Wilk y6 0.8107 <0.001 NO
Shapiro-Wilk y7 0.8715 <0.001 NO
Shapiro-Wilk y8 0.9006 <0.001 NO
Shapiro-Wilk x1 0.9751 0.1454 YES
Shapiro-Wilk x2 0.9763 0.1715 YES
Shapiro-Wilk x3 0.9734 0.1145 YES

Como era de esperarse a partir de los Q-Q plots, todas las variables \(y\) presentaron evidencia de incumplimiento del supuesto de normalidad univariada, con valores p menores a 0.05. Las variables \(x\), en cambio, no presentaron indicadores de incumplimiento de dicho supuesto.

8. Utilizar un método de estimación adecuado en caso de no cumplirse el supuesto de normalidad multivariada.

Gana y Broc ( ), establecen que para los casos de incumplimiento del supuesto de normalidad se pueden utilizar métodos de máxima verosimilitud robustos para las estimaciones, siendo el más apropiado en este caso por ser una muestra relativamente pequeña.

Otra opción es utilizar mínimos cuadrados ponderados ( ), pero no se recomienda este estimador en muestras pequeñas, sino su versión ponderada diagonalmente ( ). Este estimador, además, presenta sus versiones con errores robustos y media ajustada ( ) o con media y varianza ajustada ( ).

Por otro lado, como se mencionó en el ejercicio anterior, Forero et al. ( ) establecen que estimar con mínimos cuadrados no ponderados ( ) produce cargas factoriales más exactas y precisas que mínimos cuadrados ponderados con media y varianza ajustada ( ), por lo que sería una alternativa deseable.

El análisis de la normalidad multivariada no pareció indicar violaciones extremas al supuesto de normalidad, pero se desconoce la naturaleza de las variables originales (continuas u ordinales, cantidad de categorías, etc.). Por este motivo, la decisión más conservadora en torno a los supuestos, sería utilizar mínimos cuadrados no ponderados para la estimación del modelo.

De todas formas, también se utilizará la técnica de remuestreo (Bootstrap) para la estimación de los parámetros e indicadores de bondad de ajuste, ya que con dicha técnica no se asume ninguna distribución de los datos. La aplicación del remuestreo se empleará estimando el modelo con máxima verosimilitud robusta en 1000 iteraciones y se comparará el promedio de indicadores de bondad de ajuste y parámetros con el ajuste utilizando ULS. El modelo ajustado por máxima verosimilitud no debería considerarse por las evidencias de incumplimiento de sus supuestos.

mod <- '# Ecuaciones de medicion
             dem60 =~ y1 + y2 + y3 + y4
             dem65 =~ y5 + y6 + y7 + y8
             ind60 =~ x1 + x2 + x3
            
        # Ecuaciones estructurales
             dem60 ~ ind60
             dem65 ~ ind60 + dem60 
             
        # Correlaciones de residuos
             y1 ~~ y5
             y2 ~~ y6
             y3 ~~ y7
             y4 ~~ y8
'



# Ajuste del modelo con ULS
fit_uls <- sem(mod, data = base2, estimator="ULS")

# utilización del método de boostrap
parametros <- data.frame("dem60=~y1"=rep(NA,1000),
                         "dem60=~y2"=rep(NA,1000),
                         "dem60=~y3"=rep(NA,1000),
                         "dem60=~y4"=rep(NA,1000),
                         "dem65=~y5"=rep(NA,1000),
                         "dem65=~y6"=rep(NA,1000),
                         "dem65=~y7"=rep(NA,1000),
                         "dem65=~y8"=rep(NA,1000),
                         "ind60=~x1"=rep(NA,1000),
                         "ind60=~x2"=rep(NA,1000),
                         "ind60=~x3"=rep(NA,1000),
                         "dem60~ind60"=rep(NA,1000),
                         "dem65~ind60"=rep(NA,1000),
                         "dem65~dem60"=rep(NA,1000),
                         "y1~~y5"=rep(NA,1000),
                         "y2~~y6"=rep(NA,1000),
                         "y3~~y7"=rep(NA,1000),
                         "y4~~y8"=rep(NA,1000),
                         "y1~~y1"=rep(NA,1000),
                         "y2~~y2"=rep(NA,1000),
                         "y3~~y3"=rep(NA,1000),
                         "y4~~y4"=rep(NA,1000),
                         "y5~~y5"=rep(NA,1000),
                         "y6~~y6"=rep(NA,1000),
                         "y7~~y7"=rep(NA,1000),
                         "y8~~y8"=rep(NA,1000),
                         "x1~~x1"=rep(NA,1000),
                         "x2~~x2"=rep(NA,1000),
                         "x3~~x3"=rep(NA,1000),
                         "dem60~~dem60"=rep(NA,1000),
                         "dem65~~dem65"=rep(NA,1000),
                         "ind60~~ind60"=rep(NA,1000))
indices <- data.frame("cfi"=rep(NA,1000),
                      "tli"=rep(NA,1000),
                      "rmsea"=rep(NA,1000),
                      "srmr"=rep(NA,1000),
                      "chisq"=rep(NA,1000),
                      "df"=rep(NA,1000))
for(i in 1:1000){
  # Base de la iteracion i
  base.boot <- base2[sample(1:75,75,replace = T),]
  
  # Ajuste con la base bootstrapeada
  ajuste <- sem(mod, data = base.boot, estimator="MLR")
  
  # Guardo los parametros estandarizados en la base de parametros
  parametros[i,] <- summary(ajuste, standardized=T)$pe$std.all
  
  # Guardo los indices de bondad de ajuste
  indices$cfi[i] <- fitmeasures(ajuste)["cfi"]
  indices$tli[i] <- fitmeasures(ajuste)["tli"]
  indices$rmsea[i] <- fitmeasures(ajuste)["rmsea"]
  indices$srmr[i] <- fitmeasures(ajuste)["srmr"]
  indices$chisq[i] <- fitmeasures(ajuste)["chisq"]
  indices$df[i] <- fitmeasures(ajuste)["df"]
}




tabla_ajustes <- data.frame("Estimador"=c("ULS","Bootstrap"),
                               "CFI"=c(fitmeasures(fit_uls)[c("cfi")],
                                       mean(indices$cfi)),
                               "TLI"=c(fitmeasures(fit_uls)[c("tli")],
                                       mean(indices$tli)),
                               "RMSEA"=c(fitmeasures(fit_uls)[c("rmsea")],
                                       mean(indices$rmsea)),
                               "SRMR"=c(fitmeasures(fit_uls)[c("srmr")],
                                       mean(indices$srmr)),
                               "X2/gl"= c(fitmeasures(fit_uls)[c("chisq")] / fitmeasures(fit_uls)[c("df")],
                                          mean(indices$chisq) / mean(indices$df)),
                               check.names=F)
tabla_ajustes[,2:6] <- round(tabla_ajustes[,2:6],3)
rownames(tabla_ajustes) <- NULL

kable(tabla_ajustes,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Tabla 7. Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)
Tabla 7. Indicadores de ajuste según estimador
Estimador CFI TLI RMSEA SRMR X2/gl
ULS 0.996 0.994 0.378 0.048 11.573
Bootstrap 0.927 0.892 0.137 0.066 2.444

Ajuste con 

grafico <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,13),.85), 
          edge.label.cex = c(rep(.92,11),rep(1.2,3),rep(.87,4),rep(.75,17)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, curvature = c(1.85, 2.45, 3.05, 3.65, 1.25),
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

Ajuste con remuestreo y

Conclusiones

Habiendo descartado el ajuste con máxima verosimilitud por evidencia de incumplimiento de supuestos, se ajustó el modelo propuesto con mínimos cuadrados no ponderados (ULS) y a partir de la técnica de remuestreo tomando de base un ajuste con máxima verosimilitud robusta.

Los índices de bondad de ajuste fueron buenos en ambos escenarios (los calculados por ULS y el promedio de los 1000 ajustes por MLR con remuestreo). El promedio de los indicadores de bondad de ajuste presentó un mejor nivel en indicadores de ajuste parsimonioso (\(\Delta_{RMSEA} = 0.24\)) y en el cociente \(\chi^2 / gl\) (\(\Delta_{\chi^2/gl} = 9.106\)). El resto de índices de bondad de ajuste fueron mejores en las estimaciones con ULS (\(\Delta_{SRMR} = -0.018\), \(\Delta_{CFI} = 0.07\) y \(\Delta_{TLI} = .103\)).

Si bien no hay grandes diferencias entre un ajuste y el otro, dada la simplicidad de la estimación del modelo con ULS respecto del modelo con remuestreo y su mayor frecuencia en la literatura previa, se sugiere reportar ese modelo.

Por otro lado, los parámetros estimados por ambos caminos también son muy semejantes, por lo que se refuerza la decisión de reportar el ajuste estimado con mínimos cuadrados no ponderados (ULS)

Referencias

Forero, C. G., Maydeu-Olivares, A., & Gallardo-Pujol, D. (2009). Factor analysis with ordinal indicator: A Monte Carlo Study Comparing DWLS and ULS Estimation. Structural Equation Modeling, 16, 625–641.

Gana, K., & Broc, G. (2019). Structural equation modeling with lavaan. John Wiley & Sons.