Pasos que se van a dar.
Los datos están en un fichero dta de Stata, los podemos leer en R utilizando la librería foreign. De paso obtenemos un breve summary de las variables que nos interesan para el modelo de medida
library(foreign)
datos <- read.dta("data/paper3_seleccionvars.dta")
names(datos) <- tolower(names(datos)) # pasamos a minúsculas los nombres de las variables
Creamos vectores con los nombres de las variables que nos interesan
var.str <- c(paste("c3_", 1:3, sep = ""), "csc")
var.trust <- c(paste("d4_", 2:3, sep = ""), "d1", "d2", paste("r1_", 1:3, sep = ""))
var.confgov <- paste("d6_", 1:4, sep = "")
summary(datos[, var.str])
## c3_1 c3_2 c3_3 csc
## Min. :1.00 Min. :1.00 Min. :1.00 Min. :0.00
## 1st Qu.:1.00 1st Qu.:1.00 1st Qu.:1.00 1st Qu.:1.00
## Median :2.00 Median :1.00 Median :1.00 Median :1.00
## Mean :2.47 Mean :1.46 Mean :1.59 Mean :1.46
## 3rd Qu.:4.00 3rd Qu.:1.00 3rd Qu.:2.00 3rd Qu.:2.00
## Max. :5.00 Max. :5.00 Max. :5.00 Max. :4.00
for (i in 1:length(var.str)) hist(datos[, var.str[i]], main = var.str[i], xlab = "")
summary(datos[, var.trust])
## d4_2 d4_3 d1 d2
## Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00
## 1st Qu.:4.00 1st Qu.:3.00 1st Qu.:3.00 1st Qu.:3.00
## Median :4.00 Median :4.00 Median :4.00 Median :4.00
## Mean :4.41 Mean :3.91 Mean :3.49 Mean :3.67
## 3rd Qu.:5.00 3rd Qu.:4.00 3rd Qu.:4.00 3rd Qu.:4.00
## Max. :5.00 Max. :5.00 Max. :5.00 Max. :5.00
## r1_1 r1_2 r1_3
## Min. :1.00 Min. :1.00 Min. :1.00
## 1st Qu.:3.00 1st Qu.:3.00 1st Qu.:2.00
## Median :4.00 Median :4.00 Median :3.00
## Mean :3.69 Mean :3.68 Mean :2.88
## 3rd Qu.:4.00 3rd Qu.:4.00 3rd Qu.:4.00
## Max. :5.00 Max. :5.00 Max. :5.00
for (i in 1:length(var.trust)) hist(datos[, var.trust[i]], main = var.trust[i],
xlab = "")
summary(datos[, var.confgov])
## d6_1 d6_2 d6_3 d6_4
## Min. :1.0 Min. :1.00 Min. :1.0 Min. :1.00
## 1st Qu.:2.0 1st Qu.:2.00 1st Qu.:2.0 1st Qu.:2.00
## Median :3.0 Median :3.00 Median :3.0 Median :3.00
## Mean :3.1 Mean :2.86 Mean :3.1 Mean :2.56
## 3rd Qu.:4.0 3rd Qu.:4.00 3rd Qu.:4.0 3rd Qu.:3.00
## Max. :5.0 Max. :5.00 Max. :5.0 Max. :5.00
## NA's :287 NA's :154
for (i in 1:length(var.confgov)) hist(datos[, var.confgov[i]], main = var.confgov[i],
xlab = "")
Efectuamos un análisis factorial exploratorio sobre todas las variables, para ver si emergen los constructos de forma natural. En R hay multitud de paquetes que realizan el análisis factorial. Voy a utilizar el paquete psych.
Lo primero que hacemos es calcular la matriz de correlaciones y dibujarla ordenando sus valores
library(psych)
correlaciones <- cor(datos[, c(var.str, var.trust, var.confgov)], use = "pairwise")
cor.plot(mat.sort(correlaciones))
## Loading required package: MASS
En el gráfico se aprecia claramente el factor structural, cuyas variables presentan nula o muy poca correlación con las variables de otros factores
El factor confgov (esquina superior derecha) también se aprecia claramente y se intuyen correlaciones con el factor trust ( en el medio ), posibles regresiones a tener en cuenta.
Por último el factor trust es más difuso. Quizá sobre alguna variable o sean dos factores en vez de uno. El análisis factorial confirmatorio podría darnos más pistas
La función fa del paquete psych ofrece múltiples opciones para realizar el análisis factorial exploratorio. Veamos la más común, consistente en un análisis sin rotación mediante el método de principal axis. Elegimos 5 factores y veamos el screeplot de los valores propios para elegir un número de
factores.
fit.pa5 <- fa(datos[, c(var.str, var.trust, var.confgov)], nfactors = 5, rotate = "none",
fm = "pa")
# veamos las cargas factoriales
fit.pa5$loadings
##
## Loadings:
## PA1 PA2 PA3 PA4 PA5
## c3_1 0.168 0.540 0.183
## c3_2 0.555 0.210 0.131
## c3_3 0.483 0.285
## csc 0.666 0.318
## d4_2 0.455 0.220 -0.279 -0.183 0.428
## d4_3 0.674 0.207 -0.421 -0.273
## d1 0.495 0.118 -0.157 -0.123
## d2 0.584 0.119 -0.342 -0.209 -0.368
## r1_1 0.403 -0.186 0.363
## r1_2 0.437 0.125 -0.117 0.381
## r1_3 0.400 -0.185 0.428
## d6_1 0.587 -0.246 0.319
## d6_2 0.648 -0.234 0.398 -0.107
## d6_3 0.575 -0.332 0.321
## d6_4 0.527 -0.239 0.316
##
## PA1 PA2 PA3 PA4 PA5
## SS loadings 3.180 1.698 1.202 0.650 0.364
## Proportion Var 0.212 0.113 0.080 0.043 0.024
## Cumulative Var 0.212 0.325 0.405 0.449 0.473
plot(fit.pa5$e.values, type = "b", col = "darkred", lwd = 2, lty = 2, main = "Scree plot")
abline(h = 1)
En primer lugar vemos las cargas factoriales que nos dan una idea de en qué factor carga cada variable. También obtenemos la varianza explicada por cada factor y vemos que con 3 factores se explica el 40,5% de las variabilidad y con 4 el 44,9%.
El gráfico del análisis sería
diagram(fit.pa5, e.size = 0.1)
Se observa claramente el factor structural , pero las variables que componen los otros dos factores o constructos cargan en el mismo factor. Lo que confirma el análisis de las correlaciones.
Ahora bien, podemos utilizar rotación varimax y seleccionar 3 componentes
fit.pa4 <- fa(datos[, c(var.str, var.trust, var.confgov)], nfactors = 3, rotate = "varimax",
fm = "pa") # necesita paquete GPArotation
## Loading required package: GPArotation
diagram(fit.pa4, e.size = 0.1)
Si seleccionamos 4 factores como indicaba el scree plot y rotamos la solución tenemos
fit.pa3 <- fa(datos[, c(var.str, var.trust, var.confgov)], nfactors = 4, rotate = "varimax",
fm = "pa") # necesita paquete GPArotation
diagram(fit.pa3, e.size = 0.1)
El análisis factorial exploratorio nos indica que nuestros constructos son correctos,pero quizá haya otro más con las variables r1_1, r1_2, r1_3 ¿Tiene sentido teórico?
la función alpha del paquete psych
alpha(datos[, var.str])
##
## Reliability analysis
## Call: alpha(x = datos[, var.str])
##
## raw_alpha std.alpha G6(smc) average_r mean sd
## 0.7 0.71 0.66 0.38 1.7 0.81
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r
## c3_1 0.65 0.66 0.58 0.39
## c3_2 0.63 0.66 0.57 0.39
## c3_3 0.65 0.68 0.58 0.41
## csc 0.60 0.60 0.51 0.34
##
## Item statistics
## n r r.cor r.drop mean sd
## c3_1 998 0.72 0.57 0.48 2.5 1.38
## c3_2 998 0.73 0.58 0.49 1.5 1.00
## c3_3 998 0.70 0.55 0.44 1.6 1.13
## csc 998 0.78 0.68 0.56 1.5 0.89
##
## Non missing response frequency for each item
## 0 1 2 3 4 5 miss
## c3_1 0.00 0.34 0.24 0.14 0.19 0.10 0
## c3_2 0.00 0.78 0.09 0.04 0.06 0.03 0
## c3_3 0.00 0.73 0.10 0.05 0.07 0.04 0
## csc 0.12 0.44 0.33 0.10 0.02 0.00 0
alpha(datos[, var.trust])
##
## Reliability analysis
## Call: alpha(x = datos[, var.trust])
##
## raw_alpha std.alpha G6(smc) average_r mean sd
## 0.73 0.74 0.75 0.29 3.7 0.51
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r
## d4_2 0.72 0.72 0.71 0.30
## d4_3 0.67 0.67 0.65 0.25
## d1 0.70 0.71 0.72 0.29
## d2 0.69 0.70 0.69 0.28
## r1_1 0.71 0.72 0.72 0.30
## r1_2 0.71 0.72 0.72 0.30
## r1_3 0.72 0.73 0.73 0.31
##
## Item statistics
## n r r.cor r.drop mean sd
## d4_2 998 0.59 0.50 0.40 4.4 0.63
## d4_3 998 0.75 0.75 0.59 3.9 0.82
## d1 998 0.63 0.53 0.45 3.5 0.86
## d2 998 0.66 0.61 0.50 3.7 0.73
## r1_1 998 0.59 0.48 0.42 3.7 0.80
## r1_2 998 0.59 0.48 0.42 3.7 0.92
## r1_3 998 0.56 0.44 0.38 2.9 0.97
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## d4_2 0.00 0.00 0.06 0.46 0.48 0
## d4_3 0.00 0.04 0.25 0.46 0.25 0
## d1 0.02 0.12 0.29 0.50 0.07 0
## d2 0.00 0.07 0.25 0.60 0.08 0
## r1_1 0.01 0.11 0.15 0.65 0.08 0
## r1_2 0.01 0.14 0.13 0.58 0.13 0
## r1_3 0.03 0.41 0.24 0.29 0.03 0
alpha(datos[, var.confgov])
##
## Reliability analysis
## Call: alpha(x = datos[, var.confgov])
##
## raw_alpha std.alpha G6(smc) average_r mean sd
## 0.82 0.82 0.79 0.54 2.6 0.87
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r
## d6_1 0.78 0.78 0.71 0.55
## d6_2 0.75 0.75 0.67 0.50
## d6_3 0.77 0.77 0.71 0.53
## d6_4 0.80 0.80 0.74 0.57
##
## Item statistics
## n r r.cor r.drop mean sd
## d6_1 998 0.80 0.71 0.63 3.1 1.1
## d6_2 711 0.84 0.78 0.70 2.9 1.1
## d6_3 998 0.81 0.72 0.65 3.1 1.0
## d6_4 844 0.78 0.66 0.60 2.6 1.1
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## d6_1 0.10 0.19 0.28 0.38 0.06 0.00
## d6_2 0.14 0.23 0.30 0.31 0.03 0.29
## d6_3 0.09 0.17 0.30 0.41 0.02 0.00
## d6_4 0.19 0.29 0.29 0.21 0.01 0.15
Los \( \alpha \) son mayores que 0.7 en todas las escalas
Para el análisis factorial confirmatorio vamos a utilizar el paquete lavaan (latent variable analysis)
Primero especificamos el modelo inicial
library(lavaan)
## Loading required package: boot
## Attaching package: 'boot'
## The following object(s) are masked from 'package:psych':
##
## logit
## Loading required package: mnormt
## Loading required package: quadprog
## This is lavaan 0.5-10
## lavaan is BETA software! Please report any bugs.
# cargo algunas funciones mias útiles, que añaden por ejemplo más medidas
# de bondad de ajuste
source("../codigos_utiles/funciones_utiles.R")
modelo <- "\nstructural =~ c3_1 + c3_2 +c3_3 + csc\ntrust =~ d4_2 + d4_3 +d1 + d2 + r1_1 + r1_2 + r1_3\nconf.gov =~ d6_1 + d6_2 + d6_3 + d6_4\n"
fit.cfa1 <- cfa(modelo, data = datos, std.lv = T, std.ov = T)
# Con std.lv=T le digo que estandarice la métrica de la latente. con
# std.ov=T le digo que estandarice las variables antes del análisis es
# decir que trabaje con la matriz de correlaciones en vez de covarianzas
La especificación del modelo no se ve bien, \n indica salto de línea lo pongo tal como se ha escrito en el código
modelo <- "
structural =~ c3_1 + c3_2 +c3_3 + csc
trust =~ d4_2 + d4_3 +d1 + d2 + r1_1 + r1_2 + r1_3
conf.gov =~ d6_1 + d6_2 + d6_3 + d6_4
"
Un resumen numérico del modelo sería
summary(fit.cfa1, fit.measures = TRUE)
## lavaan (0.5-10) converged normally after 21 iterations
##
## Used Total
## Number of observations 641 998
##
## Estimator ML
## Minimum Function Chi-square 379.280
## Degrees of freedom 87
## P-value 0.000
##
## Chi-square test baseline model:
##
## Minimum Function Chi-square 2696.727
## Degrees of freedom 105
## P-value 0.000
##
## Full model versus baseline model:
##
## Comparative Fit Index (CFI) 0.887
## Tucker-Lewis Index (TLI) 0.864
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -12476.865
## Loglikelihood unrestricted model (H1) -12287.225
##
## Number of free parameters 33
## Akaike (AIC) 25019.729
## Bayesian (BIC) 25167.009
## Sample-size adjusted Bayesian (BIC) 25062.236
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.072
## 90 Percent Confidence Interval 0.065 0.080
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.056
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|)
## Latent variables:
## structural =~
## c3_1 0.558 0.043 13.087 0.000
## c3_2 0.621 0.042 14.661 0.000
## c3_3 0.566 0.043 13.285 0.000
## csc 0.766 0.042 18.100 0.000
## trust =~
## d4_2 0.532 0.040 13.274 0.000
## d4_3 0.830 0.036 22.951 0.000
## d1 0.562 0.040 14.148 0.000
## d2 0.717 0.038 19.040 0.000
## r1_1 0.395 0.042 9.505 0.000
## r1_2 0.443 0.041 10.783 0.000
## r1_3 0.326 0.042 7.744 0.000
## conf.gov =~
## d6_1 0.748 0.037 20.311 0.000
## d6_2 0.788 0.036 21.731 0.000
## d6_3 0.725 0.037 19.502 0.000
## d6_4 0.660 0.038 17.297 0.000
##
## Covariances:
## structural ~~
## trust 0.117 0.050 2.350 0.019
## conf.gov 0.040 0.050 0.789 0.430
## trust ~~
## conf.gov 0.420 0.041 10.222 0.000
##
## Variances:
## c3_1 0.687 0.046
## c3_2 0.613 0.045
## c3_3 0.678 0.046
## csc 0.411 0.046
## d4_2 0.715 0.044
## d4_3 0.310 0.033
## d1 0.683 0.042
## d2 0.485 0.036
## r1_1 0.842 0.049
## r1_2 0.802 0.047
## r1_3 0.892 0.051
## d6_1 0.439 0.034
## d6_2 0.377 0.032
## d6_3 0.473 0.034
## d6_4 0.563 0.037
## structural 1.000
## trust 1.000
## conf.gov 1.000
Vemos qeu los coeficientes sensiblemente menores de 0.6 son los asociados a las variables r1_1 r1_2 r1_3 Las medidas de ajuste serían
medidasAjuste(fit.cfa1)
## chisq df pvalue baseline.chisq
## 379.280 87.000 0.000 2696.727
## baseline.df baseline.pvalue cfi tli
## 105.000 0.000 0.887 0.864
## logl unrestricted.logl npar aic
## -12476.865 -12287.225 33.000 25019.729
## bic ntotal bic2 rmsea
## 25167.009 641.000 25062.236 0.072
## rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
## 0.065 0.080 0.000 0.056
## srmr_nomean nfi sncp ecvi
## 0.056 0.859 0.456 0.649
## gfi agfi pgfi pnfi
## 0.925 0.896 0.766 0.712
Gráficamente
# diagrama es una función mía que envuelve a diagram y a lavaan.diagram
diagrama(fit.cfa1, e.size = 0.06)
No es un modelo malo y se podría mejorar mirando los índices de modificación , que nos indican que posible relación se podría modelar para incrementar el ajuste. No obstante, vamos a probar primero sacando las variables r en otro constructo, tal como parecía indicar el factorial exploratorio y los bajos coeficientes.
modelo <- "\nstructural =~ c3_1 + c3_2 +c3_3 + csc\ntrust =~ d4_2 + d4_3 +d1 + d2\nnueva =~ r1_1 + r1_2 + r1_3\nconf.gov =~ d6_1 + d6_2 + d6_3 + d6_4\n"
fit.cfa2 <- cfa(modelo, data = datos, std.lv = T, std.ov = T)
Pongo la especificación de fit.cfa2, para que se vea mejor
modelo <- "
structural =~ c3_1 + c3_2 +c3_3 + csc
trust =~ d4_2 + d4_3 +d1 + d2
nueva =~ r1_1 + r1_2 + r1_3
conf.gov =~ d6_1 + d6_2 + d6_3 + d6_4
"
Resumen del modelo
summary(fit.cfa2, fit.measures = TRUE)
## lavaan (0.5-10) converged normally after 20 iterations
##
## Used Total
## Number of observations 641 998
##
## Estimator ML
## Minimum Function Chi-square 252.872
## Degrees of freedom 84
## P-value 0.000
##
## Chi-square test baseline model:
##
## Minimum Function Chi-square 2696.727
## Degrees of freedom 105
## P-value 0.000
##
## Full model versus baseline model:
##
## Comparative Fit Index (CFI) 0.935
## Tucker-Lewis Index (TLI) 0.919
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -12413.660
## Loglikelihood unrestricted model (H1) -12287.225
##
## Number of free parameters 36
## Akaike (AIC) 24899.321
## Bayesian (BIC) 25059.990
## Sample-size adjusted Bayesian (BIC) 24945.692
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.056
## 90 Percent Confidence Interval 0.048 0.064
## P-value RMSEA <= 0.05 0.102
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.044
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|)
## Latent variables:
## structural =~
## c3_1 0.558 0.043 13.079 0.000
## c3_2 0.620 0.042 14.643 0.000
## c3_3 0.566 0.043 13.289 0.000
## csc 0.767 0.042 18.111 0.000
## trust =~
## d4_2 0.546 0.040 13.775 0.000
## d4_3 0.881 0.036 24.434 0.000
## d1 0.534 0.040 13.430 0.000
## d2 0.711 0.038 18.816 0.000
## nueva =~
## r1_1 0.600 0.046 13.017 0.000
## r1_2 0.668 0.047 14.303 0.000
## r1_3 0.506 0.046 11.004 0.000
## conf.gov =~
## d6_1 0.747 0.037 20.287 0.000
## d6_2 0.786 0.036 21.646 0.000
## d6_3 0.727 0.037 19.580 0.000
## d6_4 0.661 0.038 17.347 0.000
##
## Covariances:
## structural ~~
## trust 0.109 0.049 2.208 0.027
## nueva 0.069 0.057 1.197 0.231
## conf.gov 0.039 0.050 0.784 0.433
## trust ~~
## nueva 0.549 0.044 12.462 0.000
## conf.gov 0.379 0.042 9.047 0.000
## nueva ~~
## conf.gov 0.372 0.050 7.464 0.000
##
## Variances:
## c3_1 0.687 0.046
## c3_2 0.614 0.045
## c3_3 0.678 0.046
## csc 0.410 0.046
## d4_2 0.701 0.043
## d4_3 0.221 0.035
## d1 0.714 0.043
## d2 0.493 0.036
## r1_1 0.638 0.050
## r1_2 0.553 0.052
## r1_3 0.742 0.050
## d6_1 0.440 0.034
## d6_2 0.381 0.032
## d6_3 0.470 0.034
## d6_4 0.561 0.037
## structural 1.000
## trust 1.000
## nueva 1.000
## conf.gov 1.000
El modelo mejora bastante, obteniendo un RMSEA con un I.C (0.048, 0.064) y valores altos de CFI y TLI Nótese que se han usado 641 observaciones (valores perdidos en d6_2) tendríamos que imputar o considerar la matriz de correlación directamente en vez de todos los datos. Queda pendiente el procedimiento de imputación, así como el modelo con las correlaciones policóricas.
Las medidas de bondad de ajuste del modelo son
medidasAjuste(fit.cfa2)
## chisq df pvalue baseline.chisq
## 252.872 84.000 0.000 2696.727
## baseline.df baseline.pvalue cfi tli
## 105.000 0.000 0.935 0.919
## logl unrestricted.logl npar aic
## -12413.660 -12287.225 36.000 24899.321
## bic ntotal bic2 rmsea
## 25059.990 641.000 24945.692 0.056
## rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
## 0.048 0.064 0.102 0.044
## srmr_nomean nfi sncp ecvi
## 0.044 0.906 0.263 0.461
## gfi agfi pgfi pnfi
## 0.950 0.928 0.760 0.725
Gráficamente
diagrama(fit.cfa2, e.size = 0.06)
Opciones que deberíamos pensar.