Cálculo de índices
Álvaro Calvo
2025-23-02
Preparación de los datos
Importamos nuestra BDD Modern Sexism
library(haven)
library(psych)## Warning: package 'psych' was built under R version 4.4.1library(ggplot2)##
## Adjuntando el paquete: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphadata <- read_sav("C:/Users/alvar/OneDrive/Escritorio/TRABAJO UDC/Analisis Estadistico Avanzado UDC/Moden Sexism.sav")Definimos nuestras variables para el cálculo de los índices, y eliminamos los valores perdidos:
myvars <- c("MS1", "MS2", "MS3", "MS4", "MS5")
data <- data[!is.na(data$MS1), ]
data <- data[!is.na(data$MS2), ]
data <- data[!is.na(data$MS3), ]
data <- data[!is.na(data$MS4), ]
data <- data[!is.na(data$MS5), ]
summary(data[, myvars])## MS1 MS2 MS3 MS4 MS5
## Min. :1.000 Min. :1.00 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:2.00 1st Qu.:1.750 1st Qu.:1.000 1st Qu.:3.000
## Median :3.000 Median :3.00 Median :3.000 Median :3.000 Median :4.000
## Mean :3.283 Mean :3.79 Mean :3.298 Mean :3.121 Mean :4.191
## 3rd Qu.:5.000 3rd Qu.:5.00 3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.:6.000
## Max. :7.000 Max. :7.00 Max. :7.000 Max. :7.000 Max. :7.000Creación de los Índices
Índice por promedio
Construimos un índice promediando la puntuación de los casos en las 5 variables de interés
indice_promedio <- rowMeans(data[, c("MS1", "MS2", "MS3", "MS4", "MS5")])
data$Indice_Promedio <- indice_promedio
summary(data$Indice_Promedio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 2.550 3.600 3.537 4.400 7.000hist(data$Indice_Promedio)
ACP Normal (psych)
suppressWarnings({library(psych)})
acp_normal <- principal(data[, c("MS1", "MS2", "MS3", "MS4", "MS5")],
nfactors = 1, rotate = "none", scores = TRUE)
puntuaciones_factoriales <- acp_normal$scores # Extraer las puntuaciones factoriales
data$PF_ACP <- puntuaciones_factoriales[, 1] # Incluir las puntuaciones factoriales del ACP a nuestra dataVemos la distribución de nuestra variable creada con un ACP normal
summary(data$PF_ACP)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.73487 -0.75133 0.01083 0.00000 0.59869 2.42122hist(data$PF_ACP)
ACP policórico (psych)
Calculamos la matriz policórica con la función polychoric del paquete psych. Recordad que usamos esta función porque nuestros datos se pueden considerar tipo ordinal. En caso de estar utilizando variables de tipo categórico, tenemos que usar la función tetrachoric de este mismo paquete (psych)
suppressWarnings({library(psych)})
data.poly <- as.data.frame(data[,c("MS1", "MS2", "MS3", "MS4", "MS5")])
r <- polychoric(data.poly)
View(r$rho)
#Creamos objeto Matriz Policórica
Matriz <- r$rhoAdecuación del modelo: KMO y Barlett
Comprobamos la adecuación del modelo factorial con el estadístico KMO y el test de esfericidad de Barlett
KMO(Matriz)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = Matriz)
## Overall MSA = 0.77
## MSA for each item =
## MS1 MS2 MS3 MS4 MS5
## 0.77 0.73 0.79 0.75 0.85cortest.bartlett(Matriz)## Warning in cortest.bartlett(Matriz): n not specified, 100 used## $chisq
## [1] 195.1124
##
## $p.value
## [1] 1.685336e-36
##
## $df
## [1] 10Decidimos cuántos factores extraer con el criterio de raíz latente
Calculamos los autovalores de cada factor que es posible extraer y se decide extraer 1 factor atendiendo al criterio de raíz latente
autovalores <- eigen(Matriz)
round(autovalores$values,2)## [1] 2.96 0.89 0.54 0.31 0.30Modelamos nuestro ACP
Especificamos la matriz de datos sobre la que calcular el modelo y el número de factores a extraer:
Modelo_poly <- pca(Matriz, nfactors = 1)Cargas factoriales y comunalidades
Vemos las cargas factoriales y las comunalidades de nuestras variables originales con el factor extraído. Ambos estadísticos son objetos ya creados dentro del Modelo.
#Observamos las cargas factoriales
print(Modelo_poly$loadings, digits = 2)##
## Loadings:
## PC1
## MS1 0.86
## MS2 0.73
## MS3 0.77
## MS4 0.80
## MS5 0.67
##
## PC1
## SS loadings 2.96
## Proportion Var 0.59#Observamos las comunalidades
print(Modelo_poly$communality, digits = 2)## MS1 MS2 MS3 MS4 MS5
## 0.74 0.53 0.60 0.64 0.46Puntuaciones factoriales
Calculamos las puntuaciones factoriales con la función factor.scores()
Puntuaciones <- factor.scores(data.poly, Modelo_poly)Para terminar, hacemos el merge de las puntuaciones calculadas con nuestra BDD original
Puntuaciones.factoriales <- Puntuaciones$scores
data <- cbind(data,Puntuaciones.factoriales)
data$PF_Poly <- data$PC1Vemos la distribución de nuestra variable creada con un ACP policórico:
summary(data$PF_Poly)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.82070 -0.77336 0.01263 0.00000 0.62273 2.53540hist(data$PF_Poly)
Comparación de modelos
Fiabilidad: Alpha de Cronbach
# Alpha para el índice por promedio
alpha_promedio <- psych::alpha(data[, c("MS1", "MS2", "MS3", "MS4", "MS5")])
print(alpha_promedio$total$raw_alpha)## [1] 0.7865588El valor alpha de Cronbach por encima de 0.70 nos señala que el índice con los promedios ya es fiable y consistente. Por tanto, seguramente las correlacioens entre el índice y las puntuaciones será elevada:
correlaciones <- cor(data[, c("Indice_Promedio", "PF_ACP", "PF_Poly")])
print(correlaciones)## Indice_Promedio PF_ACP PF_Poly
## Indice_Promedio 1.0000000 0.9992742 0.9994613
## PF_ACP 0.9992742 1.0000000 0.9998678
## PF_Poly 0.9994613 0.9998678 1.0000000Varianza explicada
acp_normal$loadings##
## Loadings:
## PC1
## MS1 0.831
## MS2 0.696
## MS3 0.739
## MS4 0.763
## MS5 0.640
##
## PC1
## SS loadings 2.713
## Proportion Var 0.543print(Modelo_poly$loadings, digits = 2)##
## Loadings:
## PC1
## MS1 0.86
## MS2 0.73
## MS3 0.77
## MS4 0.80
## MS5 0.67
##
## PC1
## SS loadings 2.96
## Proportion Var 0.59Parece que el factor policórico da cuenta de mayor varianza de las variables originales, y por tanto ajusta mejor en la reducción de la información
¿Cuál predice mejor?
#Mandamos ns y nc a perdidos
data$PTV_PP[data$PTV_PP %in% c(98, 99)] <- NA
data$PTV_VOX[data$PTV_VOX %in% c(98, 99)] <- NA
modelo_promedio <- lm(PTV_VOX ~ Indice_Promedio, data = data)
modelo_acp <- lm(PTV_VOX ~ PF_ACP, data = data)
modelo_poly <- lm(PTV_VOX ~ PF_Poly, data = data)
summary(modelo_promedio)##
## Call:
## lm(formula = PTV_VOX ~ Indice_Promedio, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2093 -1.9499 -0.7582 0.6666 9.5839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.50114 0.24937 -6.02 2.54e-09 ***
## Indice_Promedio 0.95863 0.06538 14.66 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.848 on 906 degrees of freedom
## (44 observations deleted due to missingness)
## Multiple R-squared: 0.1918, Adjusted R-squared: 0.1909
## F-statistic: 215 on 1 and 906 DF, p-value: < 2.2e-16summary(modelo_acp)##
## Call:
## lm(formula = PTV_VOX ~ PF_ACP, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2270 -1.9373 -0.7101 0.6662 9.5306
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.88859 0.09447 19.99 <2e-16 ***
## PF_ACP 1.37880 0.09375 14.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.847 on 906 degrees of freedom
## (44 observations deleted due to missingness)
## Multiple R-squared: 0.1927, Adjusted R-squared: 0.1918
## F-statistic: 216.3 on 1 and 906 DF, p-value: < 2.2e-16summary(modelo_poly)##
## Call:
## lm(formula = PTV_VOX ~ PF_Poly, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2339 -1.9417 -0.7294 0.6671 9.5271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.88899 0.09447 20.00 <2e-16 ***
## PF_Poly 1.31929 0.08968 14.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.847 on 906 degrees of freedom
## (44 observations deleted due to missingness)
## Multiple R-squared: 0.1928, Adjusted R-squared: 0.1919
## F-statistic: 216.4 on 1 and 906 DF, p-value: < 2.2e-16model4 <- lm(PTV_PP ~ Indice_Promedio, data = data)
model5 <- lm(PTV_PP ~ PF_ACP, data = data)
model6 <- lm(PTV_PP ~ PF_Poly, data = data)
summary(model4)##
## Call:
## lm(formula = PTV_PP ~ Indice_Promedio, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5625 -2.3999 -0.8186 2.1257 8.8651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02800 0.27524 0.102 0.919
## Indice_Promedio 0.79064 0.07182 11.008 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.094 on 876 degrees of freedom
## (74 observations deleted due to missingness)
## Multiple R-squared: 0.1215, Adjusted R-squared: 0.1205
## F-statistic: 121.2 on 1 and 876 DF, p-value: < 2.2e-16summary(model5)##
## Call:
## lm(formula = PTV_PP ~ PF_ACP, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5755 -2.3378 -0.8518 2.1120 8.8764
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8236 0.1044 27.05 <2e-16 ***
## PF_ACP 1.1366 0.1030 11.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.093 on 876 degrees of freedom
## (74 observations deleted due to missingness)
## Multiple R-squared: 0.122, Adjusted R-squared: 0.121
## F-statistic: 121.7 on 1 and 876 DF, p-value: < 2.2e-16summary(model6)##
## Call:
## lm(formula = PTV_PP ~ PF_Poly, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5849 -2.3566 -0.8415 2.1020 8.8732
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.82408 0.10435 27.06 <2e-16 ***
## PF_Poly 1.08892 0.09853 11.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.092 on 876 degrees of freedom
## (74 observations deleted due to missingness)
## Multiple R-squared: 0.1224, Adjusted R-squared: 0.1214
## F-statistic: 122.1 on 1 and 876 DF, p-value: < 2.2e-16Como era de esperar, el sexismo moderno predice bastante mejor el voto a VOX que al PP. En cuanto a la comparación entre índices, todos ofrecen un poder predictivo similar, si bien el policórico ajusta mínimamente mejor. En cuanto a la comparativa de los coeficientes:
coef_promedio <- data.frame(Modelo = "Promedio", Variable = names(coef(modelo_promedio)),
Coef = coef(modelo_promedio),
IC_Low = confint(modelo_promedio)[,1],
IC_High = confint(modelo_promedio)[,2])
coef_acp <- data.frame(Modelo = "ACP", Variable = names(coef(modelo_acp)),
Coef = coef(modelo_acp),
IC_Low = confint(modelo_acp)[,1],
IC_High = confint(modelo_acp)[,2])
coef_poly <- data.frame(Modelo = "Policórico", Variable = names(coef(modelo_poly)),
Coef = coef(modelo_poly),
IC_Low = confint(modelo_poly)[,1],
IC_High = confint(modelo_poly)[,2])
# Unir todo
coeficientes <- rbind(coef_promedio, coef_acp, coef_poly)
# Gráfico con intervalos de confianza
ggplot(coeficientes, aes(x = Variable, y = Coef, color = Modelo)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = IC_Low, ymax = IC_High), width = 0.2) +
labs(title = "Comparación de Coeficientes con Intervalos de Confianza",
x = "Variables",
y = "Valor del Coeficiente") +
theme_minimal()
Vemos que el ACP y el policórico tienen un coeficiente que no se diferencia.