Setting the sesion and loading the data
# Seteo de sesión
#install.packages("pacman")
library(pacman)Warning: package 'pacman' was built under R version 4.4.3pacman::p_load(
dplyr, haven, sjlabelled, psych, purrr, tidyr, sjPlot, ggplot2,
parameters, table1, car, beeswarm, lme4, scales, ggrepel, corrplot,
ggtext, patchwork, lavaan, semTools, DT, knitr, kableExtra)
options(scipen = 999)
rm(list = ls())
icils23_proc <- readRDS(here::here("~/GitHub/milenio_nudos/picils_dse/input/data/proc_data/icils23_proc.rds"))Recode the column labels
icils23_proc <- icils23_proc|>
rename(
search_info = IS3G24C,
source_info = IS3G24M,
evaluate_info = IS3G24J,
install_app = IS3G24I,
share_content = IS3G24G,
write_text = IS3G24B,
develop_webpage = IS3G24D,
create_media = IS3G24F,
insert_image = IS3G24H,
edit_image = IS3G24A,
programming = IS3G24K,
visual_coding = IS3G24L,
change_settings = IS3G24E
)The model is estimated, discarding the Safety (change_settings) and Problem solving (no items) dimensions of DigComp. Also, “Install app” was discarded because in PISA is not measured something similar (Collect and record data is in the same competence, but it is not the same task).
The CFI and TLI are excellent, indicating your model captures the latent structure well. However, the RMSEA and SRMR are high, suggesting possible model misfit or issues like item heterogeneity, correlated residuals, or threshold misspecification.
All loadings are strong (> 0.70) and highly significant (p < .001), meaning your items are good indicators of the latent constructs.
Moderate positive correlation between the two factors → related but distinct constructs.
✅ Factor structure is solid: loadings are strong, and two meaningful factors emerge. ✅ Good global fit: CFI and TLI > 0.95. ❌ Local fit concerns: RMSEA and SRMR suggest further model refinement (e.g., check residuals, thresholds, or item wording). 📊 Large sample size amplifies chi-square sensitivity—so don’t overinterpret it alone.
model_cfa <- '
gen_dse =~ search_info + source_info + evaluate_info + install_app +
share_content + write_text +
create_media + insert_image + edit_image
spec_dse =~ programming + visual_coding + develop_webpage
'
m_cfa <- cfa(model = model_cfa,
data = icils23_proc,
estimator = "DWLS",
ordered = T,
std.lv = F)
standardizedsolution(m_cfa) %>%
filter(op == "=~") %>%
select(Factor = lhs, Indicator = rhs, `Loading (DWLS)` = est.std) Factor Indicator Loading..DWLS.
1 gen_dse search_info 0.760
2 gen_dse source_info 0.670
3 gen_dse evaluate_info 0.674
4 gen_dse install_app 0.712
5 gen_dse share_content 0.747
6 gen_dse write_text 0.780
7 gen_dse create_media 0.713
8 gen_dse insert_image 0.809
9 gen_dse edit_image 0.648
10 spec_dse programming 0.819
11 spec_dse visual_coding 0.807
12 spec_dse develop_webpage 0.739summary(m_cfa)lavaan 0.6-19 ended normally after 25 iterations
Estimator DWLS
Optimization method NLMINB
Number of model parameters 49
Used Total
Number of observations 109324 135615
Model Test User Model:
Test statistic 77079.507
Degrees of freedom 53
P-value (Chi-square) 0.000
Parameter Estimates:
Parameterization Delta
Standard errors Standard
Information Expected
Information saturated (h1) model Unstructured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
gen_dse =~
search_info 1.000
source_info 0.881 0.002 393.507 0.000
evaluate_info 0.887 0.002 397.699 0.000
install_app 0.937 0.002 404.240 0.000
share_content 0.982 0.002 428.614 0.000
write_text 1.026 0.002 415.551 0.000
create_media 0.938 0.002 415.971 0.000
insert_image 1.065 0.002 437.015 0.000
edit_image 0.852 0.002 389.839 0.000
spec_dse =~
programming 1.000
visual_coding 0.985 0.003 317.303 0.000
develop_webpag 0.902 0.003 335.046 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
gen_dse ~~
spec_dse 0.389 0.001 343.981 0.000
Thresholds:
Estimate Std.Err z-value P(>|z|)
search_info|t1 -2.015 0.008 -238.223 0.000
search_info|t2 -1.434 0.006 -255.635 0.000
search_info|t3 -0.111 0.004 -29.331 0.000
source_info|t1 -1.244 0.005 -245.087 0.000
source_info|t2 -0.347 0.004 -89.606 0.000
source_info|t3 0.689 0.004 166.631 0.000
evaluate_nf|t1 -1.696 0.007 -256.299 0.000
evaluate_nf|t2 -1.022 0.005 -221.902 0.000
evaluate_nf|t3 0.369 0.004 94.986 0.000
install_app|t1 -1.838 0.007 -250.553 0.000
install_app|t2 -1.225 0.005 -243.613 0.000
install_app|t3 -0.314 0.004 -81.320 0.000
share_cntnt|t1 -1.796 0.007 -252.656 0.000
share_cntnt|t2 -1.046 0.005 -224.947 0.000
share_cntnt|t3 -0.039 0.004 -10.265 0.000
write_text|t1 -1.974 0.008 -241.549 0.000
write_text|t2 -1.379 0.005 -253.465 0.000
write_text|t3 0.014 0.004 3.817 0.000
create_medi|t1 -1.574 0.006 -257.903 0.000
create_medi|t2 -0.710 0.004 -170.782 0.000
create_medi|t3 0.369 0.004 94.974 0.000
insert_imag|t1 -1.944 0.008 -243.813 0.000
insert_imag|t2 -1.348 0.005 -251.932 0.000
insert_imag|t3 -0.257 0.004 -66.884 0.000
edit_image|t1 -1.709 0.007 -255.933 0.000
edit_image|t2 -0.953 0.004 -212.350 0.000
edit_image|t3 0.371 0.004 95.502 0.000
programming|t1 -0.868 0.004 -199.102 0.000
programming|t2 0.111 0.004 29.295 0.000
programming|t3 1.037 0.005 223.864 0.000
visual_cdng|t1 -0.874 0.004 -200.231 0.000
visual_cdng|t2 0.046 0.004 12.188 0.000
visual_cdng|t3 1.004 0.005 219.521 0.000
develp_wbpg|t1 -1.043 0.005 -224.650 0.000
develp_wbpg|t2 -0.005 0.004 -1.355 0.175
develp_wbpg|t3 0.918 0.004 207.162 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.search_info 0.422
.source_info 0.551
.evaluate_info 0.546
.install_app 0.493
.share_content 0.443
.write_text 0.392
.create_media 0.492
.insert_image 0.345
.edit_image 0.580
.programming 0.329
.visual_coding 0.349
.develop_webpag 0.454
gen_dse 0.578 0.002 311.789 0.000
spec_dse 0.671 0.003 268.435 0.000fit_pisa <- fitMeasures(m_cfa,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
knitr::kable(fit_pisa, digits = 3)| x | |
|---|---|
| chisq | 77079.50717732 |
| df | 53.00000000 |
| pvalue | 0.00000000 |
| cfi | 0.96613198 |
| tli | 0.95782473 |
| rmsea | 0.11529923 |
| srmr | 0.09026953 |
multigroup_cfa <- '
gen_dse =~ search_info + source_info + evaluate_info + install_app +
share_content + write_text +
create_media + insert_image + edit_image
spec_dse =~ programming + visual_coding + develop_webpage
'
# Configural
fitgroup <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "CNTRY",
ordered = TRUE
)
# 2 minutes
# Metrical (fix loadings)
fitgroup_metric <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "CNTRY",
ordered = TRUE,
# missing = "fiml"
group.equal = "loadings")
# 3 minutos
# Scalar (Fix loadings and intercept)
fitgroup_scalar <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "CNTRY",
ordered = TRUE,
# missing = "fiml",
group.equal = c("loadings", "intercepts"))fitmeasures_config <- fitMeasures(fitgroup,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
fitmeasures_metric <- fitMeasures(fitgroup_metric,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
fitmeasures_scalar <- fitMeasures(fitgroup_scalar,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
# Todo se demora aprox 2 minutos
# Tabla resumen (html y opción latex)
tabla_completa <- data.frame(
Modelo = c("1. Configural", "2. Métrico", "3. Escalar"),
chisq = c(fitmeasures_config["chisq"], fitmeasures_metric["chisq"], fitmeasures_scalar["chisq"]),
df = c(fitmeasures_config["df"], fitmeasures_metric["df"], fitmeasures_scalar["df"]),
CFI = c(fitmeasures_config["cfi"], fitmeasures_metric["cfi"], fitmeasures_scalar["cfi"]),
TLI = c(fitmeasures_config["tli"], fitmeasures_metric["tli"], fitmeasures_scalar["tli"]),
RMSEA = c(fitmeasures_config["rmsea"], fitmeasures_metric["rmsea"], fitmeasures_scalar["rmsea"]),
SRMR = c(fitmeasures_config["srmr"], fitmeasures_metric["srmr"], fitmeasures_scalar["srmr"])
)
# Cálculo deltas
tabla_completa$delta_chisq <- c(NA, diff(tabla_completa$chisq))
tabla_completa$delta_df <- c(NA, diff(tabla_completa$df))
tabla_completa$delta_CFI <- c(NA, diff(tabla_completa$CFI))
tabla_completa$delta_RMSEA <- c(NA, diff(tabla_completa$RMSEA))
# ANOVA para significancia
test_anova <- anova(fitgroup, fitgroup_metric, fitgroup_scalar)
# 3-4 minutos
tabla_completa$p_value_delta_chisq <- test_anova$`Pr(>Chisq)`
invariance_table_cnt <- knitr::kable(
tabla_completa,
format = "html",
digits = 3,
caption = "Análisis de invarianza por países",
col.names = c(
"Modelo",
"χ<sup>2</sup>", # chi²
"df", "CFI", "TLI", "RMSEA", "SRMR",
"Δχ<sup>2</sup>", # Δχ²
"Δdf", # Δdf
"ΔCFI", # ΔCFI
"ΔRMSEA", # ΔRMSEA
"p-valor (Δχ<sup>2</sup>)" # p-valor (Δχ²)
),
booktabs = TRUE,
escape=FALSE
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
invariance_table_cnt| Modelo | χ2 | df | CFI | TLI | RMSEA | SRMR | Δχ2 | Δdf | ΔCFI | ΔRMSEA | p-valor (Δχ2) |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1. Configural | 82176.89 | 1855 | 0.971 | 0.964 | 0.118 | 0.091 | NA | NA | NA | NA | NA |
| 2. Métrico | 89470.92 | 2195 | 0.968 | 0.967 | 0.113 | 0.094 | 7294.037 | 340 | -0.003 | -0.005 | 0 |
| 3. Escalar | 112386.68 | 2943 | 0.960 | 0.969 | 0.109 | 0.092 | 22915.759 | 748 | -0.008 | -0.004 | 0 |
icils23_proc$IS3G02 <- as.factor(icils23_proc$IS3G02)
multigroup_cfa <- '
gen_dse =~ search_info + source_info + evaluate_info + install_app +
share_content + write_text +
create_media + insert_image + edit_image
spec_dse =~ programming + visual_coding + develop_webpage
'
# Configural
fitgroup_sx <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "IS3G02",
ordered = TRUE
)Warning: lavaan->lav_data_full():
group variable 'IS3G02' contains missing values# 2 minutes
# Metrical (fix loadings)
fitgroup_metric_sx <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "IS3G02",
ordered = TRUE,
# missing = "fiml"
group.equal = "loadings")Warning: lavaan->lav_data_full():
group variable 'IS3G02' contains missing values# 3 minutos
# Scalar (Fix loadings and intercept)
fitgroup_scalar_sx <- cfa(model = multigroup_cfa,
data = icils23_proc,
group = "IS3G02",
ordered = TRUE,
# missing = "fiml",
group.equal = c("loadings", "intercepts"))Warning: lavaan->lav_data_full():
group variable 'IS3G02' contains missing values# Extracción de indicadores de ajuste
fitmeasures_config_sx <- fitMeasures(fitgroup_sx,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
fitmeasures_metric_sx <- fitMeasures(fitgroup_metric_sx,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
fitmeasures_scalar_sx <- fitMeasures(fitgroup_scalar_sx,
c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
# Tabla resumen (html y opción latex)
tabla_completa_sx <- data.frame(
Modelo = c("1. Configural", "2. Métrico", "3. Escalar"),
chisq = c(fitmeasures_config_sx["chisq"], fitmeasures_metric_sx["chisq"], fitmeasures_scalar_sx["chisq"]),
df = c(fitmeasures_config_sx["df"], fitmeasures_metric_sx["df"], fitmeasures_scalar_sx["df"]),
CFI = c(fitmeasures_config_sx["cfi"], fitmeasures_metric_sx["cfi"], fitmeasures_scalar_sx["cfi"]),
TLI = c(fitmeasures_config_sx["tli"], fitmeasures_metric_sx["tli"], fitmeasures_scalar_sx["tli"]),
RMSEA = c(fitmeasures_config_sx["rmsea"], fitmeasures_metric_sx["rmsea"], fitmeasures_scalar_sx["rmsea"]),
SRMR = c(fitmeasures_config_sx["srmr"], fitmeasures_metric_sx["srmr"], fitmeasures_scalar_sx["srmr"])
)
# Cálculo deltas
tabla_completa_sx$delta_chisq <- c(NA, diff(tabla_completa_sx$chisq))
tabla_completa_sx$delta_df <- c(NA, diff(tabla_completa_sx$df))
tabla_completa_sx$delta_CFI <- c(NA, diff(tabla_completa_sx$CFI))
tabla_completa_sx$delta_RMSEA <- c(NA, diff(tabla_completa_sx$RMSEA))
# ANOVA para significancia
test_anova_sx <- anova(fitgroup_sx, fitgroup_metric_sx, fitgroup_scalar_sx)
tabla_completa_sx$p_value_delta_chisq <- test_anova_sx$`Pr(>Chisq)`
# Tabla en formato latex
invariance_table_sx <- knitr::kable(
tabla_completa_sx,
format = "html",
digits = 3,
caption = "Análisis de invarianza por sexo",
col.names = c(
"Modelo",
"χ<sup>2</sup>", # chi²
"df", "CFI", "TLI", "RMSEA", "SRMR",
"Δχ<sup>2</sup>", # Δχ²
"Δdf", # Δdf
"ΔCFI", # ΔCFI
"ΔRMSEA", # ΔRMSEA
"p-valor (Δχ<sup>2</sup>)" # p-valor (Δχ²)
),
escape = FALSE,
booktabs = TRUE # Usar booktabs = TRUE genera tablas más profesionales en LaTeX
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
invariance_table_sx| Modelo | χ2 | df | CFI | TLI | RMSEA | SRMR | Δχ2 | Δdf | ΔCFI | ΔRMSEA | p-valor (Δχ2) |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1. Configural | 72595.35 | 106 | 0.969 | 0.962 | 0.112 | 0.088 | NA | NA | NA | NA | NA |
| 2. Métrico | 73788.89 | 116 | 0.969 | 0.964 | 0.108 | 0.088 | 1193.539 | 10 | -0.001 | -0.004 | 0 |
| 3. Escalar | 78687.72 | 138 | 0.967 | 0.968 | 0.102 | 0.088 | 4898.836 | 22 | -0.002 | -0.006 | 0 |
# Obtenemos todos los parámetros
format(Sys.time(), tz = "Chile/Continental")[1] "2025-08-06 16:26:07"all_params <- parameterEstimates(fitgroup, standardized = TRUE)
# Operador =~
factor_loadings_all_groups <- all_params %>%
filter(op == "=~")
# Extraemos las labels
group_labels <- lavInspect(fitgroup, "group.label")
# Dataframe con todos los países
country_map <- data.frame(group = 1:length(group_labels), CNT = unlist(group_labels))
# Dataframe con ambas dimensiones
factor_loadings_all_groups <- factor_loadings_all_groups %>%
left_join(country_map, by = "group")
# Dataframe gen dse
loadings_gen_dse <- factor_loadings_all_groups %>%
filter(lhs == "gen_dse") %>%
select(CNT, item = rhs, loading = std.all, latent_factor = lhs)
# Etiquetas tres países con peor ajuste
labels_gen_dse <- loadings_gen_dse %>%
group_by(item) %>%
slice_min(order_by = loading, n = 3) %>%
ungroup()
# Dataframe spec dse
loadings_spec_dse <- factor_loadings_all_groups %>%
filter(lhs == "spec_dse") %>%
select(CNT, item = rhs, loading = std.all, latent_factor = lhs)
# Etiquetas tres países con menor ajuste
labels_spec_dse <- loadings_spec_dse %>%
group_by(item) %>%
slice_min(order_by = loading, n = 3) %>%
ungroup()
# Gráfico cleaveland autoeficacia general
library(forcats)Warning: package 'forcats' was built under R version 4.4.3
Attaching package: 'forcats'The following object is masked from 'package:sjlabelled':
as_factorplot_gen_dse <- ggplot(loadings_gen_dse,
aes(x = loading,
y = fct_reorder(item, loading, .fun = median, .desc = FALSE))) +
geom_point(aes(color = '#B42012'), size = 3, alpha = 0.7) +
geom_text(data = labels_gen_dse,
aes(label = CNT),
nudge_y = 0.3,
size = 2.5,
color = "black",
check_overlap = TRUE) +
scale_color_identity(name = "Factor",
guide = "legend",
labels = c("General")) +
labs(title = "Cargas Factoriales: Autoeficacia General",
x = "Carga Factorial Estandarizada",
y = "Ítem") +
theme_minimal(base_size = 10) +
theme(axis.text.y = element_text(size = 8),
plot.title = element_text(hjust = 0.5),
legend.position = "top")
plot_gen_dse
# Cleaveland autoeficacia especializada
# Fijar el gráfico a x = .3
plot_spec_dse <- ggplot(loadings_spec_dse,
aes(x = loading,
y = fct_reorder(item, loading, .fun = median, .desc = FALSE))) +
geom_point(aes(color = '#E16462'), size = 3, alpha = 0.7) +
geom_text_repel(data = labels_spec_dse,
aes(label = CNT),
size = 2.5,
color = "black",
box.padding = 0.4,
point.padding = 0.2,
min.segment.length = 0,
max.overlaps = Inf) +
scale_color_identity(name = "Factor",
guide = "legend",
labels = c("Especializada")) +
labs(title = "Cargas Factoriales: Autoeficacia Especializada",
x = "Carga Factorial Estandarizada",
y = "Ítem") +
theme_minimal(base_size = 10) +
theme(axis.text.y = element_text(size = 8),
plot.title = element_text(hjust = 0.5),
legend.position = "top")
plot_spec_dse
# SEM por separado para cada país
lista_paises <- unique(unlist(icils23_proc$CNTRY))
resultados_lista <- list()
summary_lista <- list()
# Iteramos por cada país con el mismo modelo
for (pais_actual in lista_paises) {
format(Sys.time(), tz = "Chile/Continental")
datos_pais <- icils23_proc %>%
filter(CNTRY == pais_actual)
fit_pais <- tryCatch({
sem(model = multigroup_cfa,
data = datos_pais,
ordered = TRUE
# missing = "fiml"
)
}, error = function(e) {
cat(" -> ¡ERROR! No se pudo ajustar el modelo para", pais_actual, ". Razón:", conditionMessage(e), "\n")
return(NULL)
})
if (!is.null(fit_pais)) {
if (lavInspect(fit_pais, "converged")) {
medidas_ajuste <- fitMeasures(fit_pais)
resultados_lista[[pais_actual]] <- medidas_ajuste
} else {
cat(" -> ¡ADVERTENCIA! El modelo para", pais_actual, "no convergió. Se omitirán sus resultados.\n")
resultados_lista[[pais_actual]] <- rep(NA)
}
}
format(Sys.time(), tz = "Chile/Continental")
}
### Test con chile (o otros países)
datos_test <- icils23_proc %>%
filter(CNTRY == "CHL")
test_model <- sem(model = multigroup_cfa,
data = datos_test,
ordered = TRUE)
#fitMeasures(test_model)
# Convertimos resultados_lista en dataframe y agregamos CNT para detectar país
fit_indices_por_pais <- do.call(rbind, resultados_lista)
fit_indices_por_pais <- as.data.frame(fit_indices_por_pais)
fit_indices_por_pais <- fit_indices_por_pais %>%
tibble::rownames_to_column(var = "CNT")
# Seleccionamos indicadores de interés
indices_ajuste <- fit_indices_por_pais %>%
select(CNT, chisq, df, rmsea, rmsea.ci.lower, rmsea.ci.upper, srmr, cfi, tli)
datatable(
indices_ajuste,
options = list(
pageLength = 10,
autoWidth = TRUE,
searchHighlight = TRUE
),
filter = 'top',
rownames = FALSE,
caption = "Índices de Ajuste del Modelo (Interactiva)"
)