Economic Freedom Index
AS
2025-08-21
Setup
Libraries
Minimum set of libraries
remove(list=ls())
library(readxl)
library(tidyverse)
library(labelled)
library(ggcorrplot)
library(gridExtra)
library(GGally)
library(knitr)
library(kableExtra)
require(psych)
library(GPArotation)
require(lavaan)
require(lavaanPlot)
require(ResourceSelection)
require(semPlot)
require(tidySEM)
require(performance)
require(car)Functions
Correlation, printing, and citing
corplot <- function(d, m = "pearson", sz = 2.5, title = NULL) {
d_num <- d %>% dplyr::select(where(is.numeric))
if (ncol(d_num) < 2) return(NULL)
mycorr <- cor(d_num, use = "pairwise.complete.obs", method = m)
p.mat <- ggcorrplot::cor_pmat(d_num)
ggcorrplot(mycorr,
hc.order = TRUE,
type = "lower",
colors = c("red","white","green"),
tl.cex = 8,
tl.col = "black",
lab = TRUE,
lab_size = sz,
p.mat = p.mat,
sig.level = 0.01,
insig = "pch",
pch = 4) +
ggplot2::ggtitle(title %||% "")
}
myprint=function(x){x%>%kbl()%>%kable_classic(html_font = "Cambria")}Data
# Import Data
data <- read_excel("../data/raw/economic-freedom-of-the-world-2024-master-index-data-for-researchers-iso.xlsx",
sheet = "EFW Data 2024 Report",
range = "B5:CH4790"
)## New names:
## • `data` -> `data...9`
## • `data` -> `data...11`
## • `data` -> `data...13`
## • `data` -> `data...15`
## • `data` -> `data...17`
## • `data` -> `data...35`
## • `data` -> `data...37`
## • `data` -> `data...39`
## • `data` -> `data...44`
## • `data` -> `data...46`
## • `data` -> `data...48`df <- data
#names(df)
#str(df)Data Cleaning
Area 1: Size of Government
1A: Government consumption
1B: Transfers and subsidies
1C: Government investment
1D: Top marginal tax rate
i. 1Di: Top marginal income tax rate
ii. 1Dii: Top marginal income and payroll tax rates
1E: State ownership of assets
# --- Rename only Area 1 columns ---
df <- df %>%
rename(
`1A` = `Government consumption`, # score
`1A_data` = `data...9`, # underlying data
`1B` = `Transfers and subsidies`, # score
`1B_data` = `data...11`, # underlying data
`1C` = `Government investment`, # score
`1C_data` = `data...13`, # underlying data
`1D_i` = `Top marginal income tax rate`, # score
`1D_i_data` = `data...15`, # underlying data
`1D_ii` = `Top marginal income and payroll tax rate`, # score
`1D_ii_data` = `data...17`, # underlying data
`1D` = `Top marginal tax rate`, # combined top rate
`1E` = `State ownership of Assets`, # score
`Area1` = `Size of Government`, # Area 1 summary score
`Area1_rank` = `Area 1 Rank` # Area 1 rank
)
# --- Attach labels so you keep the readable names ---
var_label(df$`1A`) <- "Government consumption"
var_label(df$`1A_data`) <- "Government consumption – data"
var_label(df$`1B`) <- "Transfers and subsidies"
var_label(df$`1B_data`) <- "Transfers and subsidies – data"
var_label(df$`1C`) <- "Government investment"
var_label(df$`1C_data`) <- "Government investment – data"
var_label(df$`1D_i`) <- "Top marginal income tax rate"
var_label(df$`1D_i_data`) <- "Top marginal income tax rate – data"
var_label(df$`1D_ii`) <- "Top marginal income and payroll tax rate"
var_label(df$`1D_ii_data`) <- "Top marginal income and payroll tax rate – data"
var_label(df$`1D`) <- "Top marginal tax rate (aggregate score)"
var_label(df$`1E`) <- "State ownership of Assets"
var_label(df$Area1) <- "Area 1: Size of Government (summary score)"
var_label(df$Area1_rank) <- "Area 1 Rank"
#dplyr::glimpse(df) # see new names
#sapply(df[c("1A","1B","1C","1D","1E","Area1")], var_label) # see labelsArea 2: Legal System and Property Rights
2A–2H cover the eight sub-indicators.
2_adj, 2_wo_gen, 2_with_gen capture the gender adjustment and the two alternate summary indices.
Area2 = the official “with gender adjustment” score.
df <- df %>%
rename(
`2A` = `Judicial independence`,
`2B` = `Impartial courts`,
`2C` = `Protection of property rights`,
`2D` = `Military interference in rule of law and politics`,
`2E` = `Integrity of the legal system`,
`2F` = `Legal enforcement of contracts`,
`2G` = `Regulatory restrictions on the sale of real property`,
`2H` = `Reliability of police`,
`2_adj` = `Gender Legal Rights Adjustment`,
`2_wo_gen` = `Legal System & Property Rights -- Without Gender Adjustment`,
`2_with_gen` = `Legal System & Property Rights -- With Gender Adjustment`,
Area2_rank = `Area 2 Rank`
) %>%
mutate(
Area2 = `2_with_gen`
) %>%
relocate(Area2, .after = `2_with_gen`)
# Labels
var_label(df$`2A`) <- "Judicial independence"
var_label(df$`2B`) <- "Impartial courts"
var_label(df$`2C`) <- "Protection of property rights"
var_label(df$`2D`) <- "Military interference in rule of law and politics"
var_label(df$`2E`) <- "Integrity of the legal system"
var_label(df$`2F`) <- "Legal enforcement of contracts"
var_label(df$`2G`) <- "Regulatory restrictions on the sale of real property"
var_label(df$`2H`) <- "Reliability of police"
var_label(df$`2_adj`) <- "Gender Legal Rights Adjustment"
var_label(df$`2_wo_gen`) <- "Legal System & Property Rights (without Gender Adjustment)"
var_label(df$`2_with_gen`) <- "Legal System & Property Rights (with Gender Adjustment)"
var_label(df$Area2) <- "Area 2: Legal System & Property Rights (summary score)"
var_label(df$Area2_rank) <- "Area 2: Rank"Area 3: Sound Money
df <- df %>%
rename(
`3A` = `Money growth`,
`3A_data` = `data...35`,
`3B` = `Standard deviation of inflation`,
`3B_data` = `data...37`,
`3C` = `Inflation: Most recent year`,
`3C_data` = `data...39`,
`3D` = `Freedom to own foreign currency bank accounts`,
Area3 = `Sound Money`,
Area3_rank = `Area 3 Rank`
)
var_label(df$`3A`) <- "Money growth"
var_label(df$`3A_data`) <- "Money growth – data"
var_label(df$`3B`) <- "Standard deviation of inflation"
var_label(df$`3B_data`) <- "Standard deviation of inflation – data"
var_label(df$`3C`) <- "Inflation: Most recent year"
var_label(df$`3C_data`) <- "Inflation: Most recent year – data"
var_label(df$`3D`) <- "Freedom to own foreign currency bank accounts"
var_label(df$Area3) <- "Area 3: Sound Money (summary score)"
var_label(df$Area3_rank) <- "Area 3: Rank"Area 4: Freedom to Trade Internationally
df <- df %>%
rename(
`4A_i` = `Revenue from trade taxes (% of trade sector)`,
`4A_i_data` = `data...44`,
`4A_ii` = `Mean tariff rate`,
`4A_ii_data`= `data...46`,
`4A_iii` = `Standard deviation of tariff rates`,
`4A_iii_data`= `data...48`,
`4A` = `Tariffs`,
`4B_i` = `Non-tariff trade barriers`,
`4B_ii` = `Compliance costs of importing and exporting`,
`4B` = `Regulatory trade barriers`,
`4C` = `Black market exchange rates`,
`4D_i` = `Financial openness`,
`4D_ii` = `Capital controls`,
`4D_iii` = `Freedom of foreigners to visit`,
`4D_iv` = `Protection of foreign assets`,
`4D` = `Controls of the movement of capital and people`,
Area4 = `Freedom to trade internationally`,
Area4_rank = `Area 4 Rank`
)
var_label(df$`4A_i`) <- "Revenue from trade taxes (% of trade sector)"
var_label(df$`4A_i_data`) <- "Revenue from trade taxes – data"
var_label(df$`4A_ii`) <- "Mean tariff rate"
var_label(df$`4A_ii_data`)<- "Mean tariff rate – data"
var_label(df$`4A_iii`) <- "Standard deviation of tariff rates"
var_label(df$`4A_iii_data`)<- "Standard deviation of tariff rates – data"
var_label(df$`4A`) <- "Tariffs (aggregate score)"
var_label(df$`4B_i`) <- "Non-tariff trade barriers"
var_label(df$`4B_ii`) <- "Compliance costs of importing and exporting"
var_label(df$`4B`) <- "Regulatory trade barriers (aggregate score)"
var_label(df$`4C`) <- "Black market exchange rates"
var_label(df$`4D_i`) <- "Financial openness"
var_label(df$`4D_ii`) <- "Capital controls"
var_label(df$`4D_iii`) <- "Freedom of foreigners to visit"
var_label(df$`4D_iv`) <- "Protection of foreign assets"
var_label(df$`4D`) <- "Controls of the movement of capital and people (aggregate score)"
var_label(df$Area4) <- "Area 4: Freedom to trade internationally (summary score)"
var_label(df$Area4_rank) <- "Area 4: Rank"Area 5: Regulation
df <- df %>%
rename(
`5A_i` = `Ownership of banks`,
`5A_ii` = `Private sector credit`,
`5A_iii` = `Interest rate controls/negative real interest rates)`,
`5A` = `Credit market regulations`,
`5B_i` = `Hiring regulations and minimum wage`,
`5B_ii` = `Hiring and firing regulations`,
`5B_iii` = `Centralized collective bargaining`,
`5B_iv` = `Hours Regulations`,
`5B_v` = `Mandated cost of worker dismissal`,
`5B_vi` = `Conscription`,
`5B_vii` = `Foreign Labor`,
`5B` = `Labor market regulations`,
`5C_i` = `Regulatory Burden`,
`5C_ii` = `Bureacracy costs`,
`5C_iii` = `Impartial Public Administration`,
`5C_iv` = `Tax compliance`,
`5C` = `Business regulations`,
`5D_i` = `Market openness`,
`5D_ii` = `Business Permits`,
`5D_iii` = `Distorton of the business environment`,
`5D` = `Freedom to enter markets and compete`,
Area5 = `Regulation`,
Area5_rank = `Area 5 Rank`
)
var_label(df$`5A_i`) <- "Ownership of banks"
var_label(df$`5A_ii`) <- "Private sector credit"
var_label(df$`5A_iii`) <- "Interest rate controls / negative real interest rates"
var_label(df$`5A`) <- "Credit market regulations (aggregate score)"
var_label(df$`5B_i`) <- "Hiring regulations and minimum wage"
var_label(df$`5B_ii`) <- "Hiring and firing regulations"
var_label(df$`5B_iii`) <- "Centralized collective bargaining"
var_label(df$`5B_iv`) <- "Hours Regulations"
var_label(df$`5B_v`) <- "Mandated cost of worker dismissal"
var_label(df$`5B_vi`) <- "Conscription"
var_label(df$`5B_vii`) <- "Foreign Labor"
var_label(df$`5B`) <- "Labor market regulations (aggregate score)"
var_label(df$`5C_i`) <- "Regulatory Burden"
var_label(df$`5C_ii`) <- "Bureaucracy costs"
var_label(df$`5C_iii`) <- "Impartial Public Administration"
var_label(df$`5C_iv`) <- "Tax compliance"
var_label(df$`5C`) <- "Business regulations (aggregate score)"
var_label(df$`5D_i`) <- "Market openness"
var_label(df$`5D_ii`) <- "Business Permits"
var_label(df$`5D_iii`) <- "Distortion of the business environment"
var_label(df$`5D`) <- "Freedom to enter markets and compete (aggregate score)"
var_label(df$Area5) <- "Area 5: Regulation (summary score)"
var_label(df$Area5_rank) <- "Area 5: Rank"library(writexl)
# export
write_xlsx(x = df, path = "economic_freedom_cleaned.xlsx")Descriptive
Descriptive prior to normalizing.
mydata <- df %>%
filter(Year == 2022) %>%
select(
Year, Countries, `Economic Freedom Summary Index`, Rank,
# Area 1
`1A`, `1B`, `1C`, `1D_i`, `1D_ii`, `1D`, `1E`, Area1, Area1_rank,
# Area 2
`2A`, `2B`, `2C`, `2D`, `2E`, `2F`, `2G`, `2H`, `2_adj`, `2_wo_gen`, `2_with_gen`, Area2, Area2_rank,
# Area 3
`3A`, `3B`, `3C`, `3D`, Area3, Area3_rank,
# Area 4
`4A_i`, `4A_ii`, `4A_iii`, `4A`, `4B_i`, `4B_ii`, `4B`, `4C`,
`4D_i`, `4D_ii`, `4D_iii`, `4D_iv`, `4D`, Area4, Area4_rank,
# Area 5
`5A_i`, `5A_ii`, `5A_iii`, `5A`,
`5B_i`, `5B_ii`, `5B_iii`, `5B_iv`, `5B_v`, `5B_vi`, `5B_vii`, `5B`,
`5C_i`, `5C_ii`, `5C_iii`, `5C_iv`, `5C`,
`5D_i`, `5D_ii`, `5D_iii`, `5D`,
Area5, Area5_rank
)
# myprint(describe(mydata))Clean Summary Stats Table
# --- Build "code: label" vector with minor corrections ---
lab_vec <- vapply(
names(mydata),
function(n) {
lb <- attr(mydata[[n]], "label", exact = TRUE)
if (is.null(lb) || is.na(lb) || lb == "") {
return(n)
}
# Clean up Area summaries and ranks
if (grepl("^Area[0-9]+$", n)) {
lb <- sub("^Area [0-9]+: ?", "", lb)
}
if (grepl("^Area[0-9]+_rank$", n)) {
lb <- "Rank"
}
paste0(n, ": ", lb)
},
character(1)
)
# --- Run describe on numeric vars ---
desc <- psych::describe(select(mydata, where(is.numeric)), fast = FALSE)
# --- Relabel rownames with code:label ---
rownames(desc) <- lab_vec[rownames(desc)]
# --- Convert to data.frame and select key stats ---
desc_df <- desc %>%
as.data.frame() %>%
rownames_to_column("Variable") %>%
select(Variable, n, mean, sd, min, max ) # median, skew, kurtosis, se
# --- Pretty table output ---
kable(
desc_df,
digits = 3,
col.names = c("Variable", "N", "Mean", "SD",
"Min", "Max"), # "Median","Skew", "Kurtosis", "SE"
align = "lccccc",
caption = "Summary Statistics with Labels"
) %>%
kable_styling(full_width = FALSE, position = "center", font_size = 12)| Variable | N | Mean | SD | Min | Max |
|---|---|---|---|---|---|
| Year | 165 | 2022.000 | 0.000 | 2022.000 | 2022.000 |
| Economic Freedom Summary Index | 165 | 6.529 | 1.049 | 3.020 | 8.580 |
| Rank | 165 | 82.776 | 47.801 | 1.000 | 165.000 |
| 1A: Government consumption | 165 | 5.471 | 2.456 | 0.000 | 10.000 |
| 1B: Transfers and subsidies | 156 | 7.526 | 2.054 | 1.823 | 10.000 |
| 1C: Government investment | 161 | 7.119 | 3.275 | 0.000 | 10.000 |
| 1D_i: Top marginal income tax rate | 165 | 7.473 | 2.194 | 2.000 | 10.000 |
| 1D_ii: Top marginal income and payroll tax rate | 162 | 5.481 | 2.633 | 0.000 | 10.000 |
| 1D: Top marginal tax rate (aggregate score) | 165 | 6.497 | 2.227 | 1.000 | 10.000 |
| 1E: State ownership of Assets | 162 | 6.664 | 1.580 | 2.435 | 9.520 |
| Area1: Size of Government (summary score) | 165 | 6.640 | 1.121 | 3.622 | 9.061 |
| Area1_rank: Rank | 165 | 83.000 | 47.776 | 1.000 | 165.000 |
| 2A: Judicial independence | 165 | 5.514 | 1.447 | 2.450 | 8.636 |
| 2B: Impartial courts | 165 | 4.660 | 1.888 | 0.827 | 8.844 |
| 2C: Protection of property rights | 165 | 5.437 | 2.157 | 0.000 | 9.668 |
| 2D: Military interference in rule of law and politics | 138 | 6.292 | 2.660 | 0.000 | 10.000 |
| 2E: Integrity of the legal system | 164 | 5.830 | 1.848 | 1.704 | 9.797 |
| 2F: Legal enforcement of contracts | 165 | 3.975 | 1.983 | 0.000 | 8.728 |
| 2G: Regulatory restrictions on the sale of real property | 163 | 7.591 | 1.522 | 2.666 | 9.981 |
| 2H: Reliability of police | 165 | 5.206 | 2.420 | 0.000 | 9.769 |
| 2_adj: Gender Legal Rights Adjustment | 165 | 0.871 | 0.167 | 0.294 | 1.000 |
| 2_wo_gen: Legal System & Property Rights (without Gender Adjustment) | 165 | 5.525 | 1.660 | 1.984 | 9.101 |
| 2_with_gen: Legal System & Property Rights (with Gender Adjustment) | 165 | 5.227 | 1.790 | 1.634 | 9.101 |
| Area2: Legal System & Property Rights (summary score) | 165 | 5.227 | 1.790 | 1.634 | 9.101 |
| Area2_rank: Rank | 165 | 83.000 | 47.776 | 1.000 | 165.000 |
| 3A: Money growth | 163 | 8.242 | 1.614 | 0.000 | 9.979 |
| 3B: Standard deviation of inflation | 165 | 7.652 | 2.520 | 0.000 | 9.740 |
| 3C: Inflation: Most recent year | 165 | 5.885 | 2.610 | 0.000 | 9.460 |
| 3D: Freedom to own foreign currency bank accounts | 165 | 7.152 | 3.952 | 0.000 | 10.000 |
| Area3: Sound Money (summary score) | 165 | 7.235 | 1.712 | 0.739 | 9.552 |
| Area3_rank: Rank | 165 | 82.994 | 47.765 | 1.000 | 165.000 |
| 4A_i: Revenue from trade taxes (% of trade sector) | 156 | 8.631 | 1.469 | 3.327 | 10.000 |
| 4A_ii: Mean tariff rate | 164 | 8.261 | 1.017 | 3.500 | 10.000 |
| 4A_iii: Standard deviation of tariff rates | 165 | 6.057 | 2.127 | 0.000 | 10.000 |
| 4A: Tariffs (aggregate score) | 165 | 7.621 | 1.177 | 3.917 | 10.000 |
| 4B_i: Non-tariff trade barriers | 165 | 5.754 | 1.804 | 0.000 | 9.182 |
| 4B_ii: Compliance costs of importing and exporting | 165 | 6.456 | 3.010 | 0.000 | 9.982 |
| 4B: Regulatory trade barriers (aggregate score) | 165 | 6.105 | 2.112 | 1.688 | 9.414 |
| 4C: Black market exchange rates | 165 | 9.219 | 2.474 | 0.000 | 10.000 |
| 4D_i: Financial openness | 165 | 5.844 | 3.016 | 0.000 | 10.000 |
| 4D_ii: Capital controls | 164 | 3.516 | 2.758 | 0.000 | 10.000 |
| 4D_iii: Freedom of foreigners to visit | 165 | 7.104 | 3.380 | 0.000 | 10.000 |
| 4D_iv: Protection of foreign assets | 165 | 5.727 | 2.099 | 0.000 | 9.266 |
| 4D: Controls of the movement of capital and people (aggregate score) | 165 | 5.553 | 2.145 | 0.690 | 9.400 |
| Area4: Freedom to trade internationally (summary score) | 165 | 7.125 | 1.490 | 2.480 | 9.661 |
| Area4_rank: Rank | 165 | 83.000 | 47.776 | 1.000 | 165.000 |
| 5A_i: Ownership of banks | 151 | 7.589 | 2.784 | 0.000 | 10.000 |
| 5A_ii: Private sector credit | 162 | 7.392 | 2.535 | 0.000 | 10.000 |
| 5A_iii: Interest rate controls / negative real interest rates | 163 | 8.350 | 2.268 | 0.000 | 10.000 |
| 5A: Credit market regulations (aggregate score) | 164 | 7.743 | 1.631 | 0.000 | 10.000 |
| 5B_i: Hiring regulations and minimum wage | 165 | 5.815 | 1.800 | 1.250 | 9.450 |
| 5B_ii: Hiring and firing regulations | 165 | 5.642 | 1.969 | 0.000 | 10.000 |
| 5B_iii: Centralized collective bargaining | 141 | 6.481 | 1.220 | 1.611 | 8.672 |
| 5B_iv: Hours Regulations | 165 | 7.842 | 1.938 | 2.000 | 10.000 |
| 5B_v: Mandated cost of worker dismissal | 162 | 6.494 | 2.835 | 0.000 | 10.000 |
| 5B_vi: Conscription | 165 | 6.612 | 4.203 | 0.000 | 10.000 |
| 5B_vii: Foreign Labor | 165 | 5.018 | 1.495 | 0.000 | 8.139 |
| 5B: Labor market regulations (aggregate score) | 165 | 6.254 | 1.221 | 2.674 | 9.137 |
| 5C_i: Regulatory Burden | 158 | 4.148 | 1.223 | 1.137 | 7.439 |
| 5C_ii: Bureaucracy costs | 165 | 5.036 | 2.216 | 0.000 | 9.556 |
| 5C_iii: Impartial Public Administration | 162 | 5.617 | 2.482 | 0.338 | 9.910 |
| 5C_iv: Tax compliance | 165 | 5.901 | 1.924 | 0.000 | 9.874 |
| 5C: Business regulations (aggregate score) | 165 | 5.168 | 1.549 | 0.597 | 9.051 |
| 5D_i: Market openness | 165 | 6.070 | 1.795 | 1.243 | 9.995 |
| 5D_ii: Business Permits | 161 | 8.386 | 1.141 | 4.978 | 9.998 |
| 5D_iii: Distortion of the business environment | 165 | 5.201 | 2.290 | 0.000 | 10.000 |
| 5D: Freedom to enter markets and compete (aggregate score) | 165 | 6.515 | 1.436 | 2.592 | 9.856 |
| Area5: Regulation (summary score) | 165 | 6.415 | 1.127 | 2.541 | 8.855 |
| Area5_rank: Rank | 165 | 83.000 | 47.776 | 1.000 | 165.000 |
Correlation plots by Area
# ---- Helper that ignores missing vars gracefully ----
sel_vars <- function(df, vars) dplyr::select(df, dplyr::any_of(vars))
# ============
# GROUPS by Area (matches your current columns)
# ============
groups <- list(
A1 = c("1A","1B","1C","1D_i","1D_ii","1D","1E","Area1"),
A2 = c("2A","2B","2C","2D","2E","2F","2G","2H","2_adj","2_wo_gen","2_with_gen","Area2"),
A3 = c("3A","3B","3C","3D","Area3"),
A4 = c("4A_i","4A_ii","4A_iii","4A","4B_i","4B_ii","4B","4C","4D_i","4D_ii","4D_iii","4D_iv","4D","Area4"),
A5 = c("5A_i","5A_ii","5A_iii","5A",
"5B_i","5B_ii","5B_iii","5B_iv","5B_v","5B_vi","5B_vii","5B",
"5C_i","5C_ii","5C_iii","5C_iv","5C",
"5D_i","5D_ii","5D_iii","5D","Area5")
)
# Optional: keep rank columns for reference or plotting separately
ranks <- c("Area1_rank","Area2_rank","Area3_rank","Area4_rank","Area5_rank")
# ============
# EXAMPLE: build a 2x3 grid of correlation plots by Area
# (skips a plot if not enough numeric columns)
# ============
plots <- list(
corplot(sel_vars(mydata, groups$A1), title = "Area 1 (1A–1E + Area1)"),
corplot(sel_vars(mydata, groups$A2), title = "Area 2 (2A–2H + adj & summaries)"),
corplot(sel_vars(mydata, groups$A3), title = "Area 3 (3A–3D + Area3)"),
corplot(sel_vars(mydata, groups$A4), title = "Area 4 (4A–4D + Area4)"),
corplot(sel_vars(mydata, groups$A5), title = "Area 5 (5A–5D + Area5)")
)
# Drop NULLs (in case a group had <2 numeric vars)
plots <- Filter(Negate(is.null), plots)
# Arrange (2 columns looks small)
# if (length(plots) > 0) gridExtra::grid.arrange(grobs = plots, ncol = 1)
# Print each correlation plot separately
for (i in seq_along(plots)) {
print(plots[[i]])
}
all_vars <- unique(unlist(groups))
plot_all <- corplot(sel_vars(mydata, all_vars), sz = 2, title = "All Areas: indicators + Area summaries")
if (!is.null(plot_all)) print(plot_all)
KDE Pairs Plot 1
KDE Pairs of data
#dev.new(width = 12, height = 12)
# robust KDE/pairs that avoids density; uses hex -> points
kdepairs <- function(d, title = NULL) {
d_num <- d %>%
dplyr::select(where(is.numeric)) %>%
mutate(across(everything(), ~ ifelse(is.finite(.), ., NA_real_))) %>% # Inf -> NA
tidyr::drop_na()
# drop constant columns
keep <- vapply(d_num, function(x) {x <- stats::na.omit(x); length(unique(x)) > 1}, logical(1))
d_num <- d_num[, keep, drop = FALSE]
if (ncol(d_num) < 2L || nrow(d_num) < 3L) return(NULL)
make_plot <- function(lower_panel) {
GGally::ggpairs(
d_num,
lower = list(continuous = lower_panel),
diag = list(continuous = GGally::wrap("densityDiag")),
upper = list(continuous = GGally::wrap("cor", size = 3))
) +
labs(title = title) +
theme(strip.text = element_text(size = 9),
axis.text.x = element_text(size = 8, angle = 45, hjust = 1),
axis.text.y = element_text(size = 8))
}
# Use hex by default (robust to NA/Inf after cleaning)
p <- tryCatch(make_plot(GGally::wrap("hex", bins = 20)),
error = function(e) NULL)
if (is.null(p)) {
# Final fallback: points (never calls kde2d / density2d)
p <- make_plot(GGally::wrap("points", alpha = 0.5))
}
p
}
for (grp in names(groups)) {
d_sub <- sel_vars(mydata, groups[[grp]])
p <- kdepairs(d_sub, title = paste(grp, "block"))
if (!is.null(p)) {
cat("\n\n###", grp, "block\n")
print(p)
} else {
message(sprintf("Skipping %s: fewer than 2 numeric columns after cleaning.", grp))
}
}##
##
## ### A1 block
##
##
## ### A2 block
##
##
## ### A3 block
##
##
## ### A4 block
##
##
## ### A5 block
Normalize
Z-scores to adjust for magnitudes
# 1) drop overall index/ranks and Area summaries/ranks
mydata_clean <- mydata %>%
select(
-`Economic Freedom Summary Index`,
-Rank,
-matches("^Area[1-5]$"),
-matches("^Area[1-5]_rank$")
)
# 2) standardize all numeric columns except Year (Countries is chr, so untouched)
num_cols <- mydata_clean %>%
select(where(is.numeric)) %>%
names() %>%
setdiff("Year")
mydata_std <- mydata_clean %>%
mutate(across(all_of(num_cols), ~ as.numeric(scale(.))))mydata_std has Year, Countries, all other numeric
columns z-scored, labels preserved
Descriptive of Standardized Data
myprint(describe(mydata_std))| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Year | 1 | 165 | 2022 | 0.00000 | 2022.0000000 | 2022.0000000 | 0.0000000 | 2022.000000 | 2022.0000000 | 0.000000 | NaN | NaN | 0.0000000 |
| Countries | 2 | 165 | 83 | 47.77552 | 83.0000000 | 83.0000000 | 60.7866000 | 1.000000 | 165.0000000 | 164.000000 | 0.0000000 | -1.2218392 | 3.7193189 |
| 1A | 3 | 165 | 0 | 1.00000 | 0.0523277 | 0.0443824 | 1.0989245 | -2.227595 | 1.8436956 | 4.071291 | -0.3235424 | -0.5053765 | 0.0778499 |
| 1B | 4 | 156 | 0 | 1.00000 | 0.3288157 | 0.1063876 | 0.9323013 | -2.776855 | 1.2043998 | 3.981254 | -0.8131572 | -0.2875138 | 0.0800641 |
| 1C | 5 | 161 | 0 | 1.00000 | 0.3719076 | 0.1581273 | 0.7527236 | -2.173597 | 0.8796127 | 3.053209 | -1.0410300 | -0.2304465 | 0.0788110 |
| 1D_i | 6 | 165 | 0 | 1.00000 | 0.2403719 | 0.0826999 | 1.3517687 | -2.494895 | 1.1521274 | 3.647022 | -0.4837939 | -0.8679662 | 0.0778499 |
| 1D_ii | 7 | 162 | 0 | 1.00000 | -0.1828412 | -0.0104944 | 1.1260259 | -2.081576 | 1.7158940 | 3.797470 | 0.0551937 | -0.8775455 | 0.0785674 |
| 1D | 8 | 165 | 0 | 1.00000 | 0.0013608 | 0.0233076 | 0.9986781 | -2.468503 | 1.5730922 | 4.041595 | -0.2289941 | -0.7543954 | 0.0778499 |
| 1E | 9 | 162 | 0 | 1.00000 | 0.2963537 | 0.0924341 | 0.7296157 | -2.676932 | 1.8075232 | 4.484455 | -0.8308829 | -0.0998867 | 0.0785674 |
| 2A | 10 | 165 | 0 | 1.00000 | -0.0955808 | -0.0339725 | 1.0484079 | -2.117834 | 2.1581993 | 4.276033 | 0.2833011 | -0.6546634 | 0.0778499 |
| 2B | 11 | 165 | 0 | 1.00000 | -0.0878123 | -0.0536834 | 1.0665836 | -2.029734 | 2.2156250 | 4.245359 | 0.4152645 | -0.5649875 | 0.0778499 |
| 2C | 12 | 165 | 0 | 1.00000 | -0.1413866 | -0.0179148 | 1.0488531 | -2.520281 | 1.9616347 | 4.481916 | 0.1927469 | -0.6560510 | 0.0778499 |
| 2D | 13 | 138 | 0 | 1.00000 | -0.1725161 | 0.0316577 | 1.3932860 | -2.365286 | 1.3937482 | 3.759034 | -0.1117562 | -1.1084036 | 0.0851257 |
| 2E | 14 | 164 | 0 | 1.00000 | 0.0516941 | 0.0054411 | 1.0471520 | -2.233394 | 2.1471235 | 4.380518 | -0.0188837 | -0.5990356 | 0.0780869 |
| 2F | 15 | 165 | 0 | 1.00000 | -0.1293499 | -0.0447961 | 0.9736318 | -2.004666 | 2.3974799 | 4.402146 | 0.3330804 | -0.5000382 | 0.0778499 |
| 2G | 16 | 163 | 0 | 1.00000 | 0.0563480 | 0.0813794 | 0.8535221 | -3.236745 | 1.5706358 | 4.807380 | -0.7517940 | 0.4735554 | 0.0783260 |
| 2H | 17 | 165 | 0 | 1.00000 | -0.0490894 | 0.0030132 | 1.1697726 | -2.151028 | 1.8858288 | 4.036856 | -0.0369216 | -0.8335076 | 0.0778499 |
| 2_adj | 18 | 165 | 0 | 1.00000 | 0.4205545 | 0.2007349 | 0.5222327 | -3.454098 | 0.7727957 | 4.226893 | -1.4906346 | 1.3361288 | 0.0778499 |
| 2_wo_gen | 19 | 165 | 0 | 1.00000 | -0.1946054 | -0.0332058 | 1.0155411 | -2.133918 | 2.1545685 | 4.288487 | 0.2642647 | -0.6606236 | 0.0778499 |
| 2_with_gen | 20 | 165 | 0 | 1.00000 | -0.1322642 | -0.0322780 | 1.0147784 | -2.007770 | 2.1643176 | 4.172087 | 0.2917051 | -0.6190112 | 0.0778499 |
| 3A | 21 | 163 | 0 | 1.00000 | 0.1835545 | 0.1664816 | 0.5984660 | -5.107349 | 1.0766004 | 6.183949 | -3.0279091 | 11.6825191 | 0.0783260 |
| 3B | 22 | 165 | 0 | 1.00000 | 0.4023375 | 0.2209118 | 0.4280087 | -3.037081 | 0.8286547 | 3.865735 | -1.8813173 | 2.6762947 | 0.0778499 |
| 3C | 23 | 165 | 0 | 1.00000 | 0.3093835 | 0.1308643 | 0.7021703 | -2.254840 | 1.3700188 | 3.624859 | -1.0752626 | 0.1576541 | 0.0778499 |
| 3D | 24 | 165 | 0 | 1.00000 | 0.7207972 | 0.1309911 | 0.0000000 | -1.809661 | 0.7207972 | 2.530458 | -0.9132408 | -0.8011341 | 0.0778499 |
| 4A_i | 25 | 156 | 0 | 1.00000 | 0.3831277 | 0.1860429 | 0.4854727 | -3.611203 | 0.9323482 | 4.543552 | -1.7252757 | 2.6904185 | 0.0800641 |
| 4A_ii | 26 | 164 | 0 | 1.00000 | 0.2646989 | 0.1043625 | 0.8165365 | -4.682183 | 1.7104080 | 6.392591 | -1.1335645 | 2.0006803 | 0.0780869 |
| 4A_iii | 27 | 165 | 0 | 1.00000 | 0.3671546 | 0.1266397 | 0.4925761 | -2.847112 | 1.8534828 | 4.700595 | -1.2978652 | 1.7794047 | 0.0778499 |
| 4A | 28 | 165 | 0 | 1.00000 | 0.2292098 | 0.0703047 | 0.7646174 | -3.148247 | 2.0219768 | 5.170224 | -0.6675823 | 0.0321753 | 0.0778499 |
| 4B_i | 29 | 165 | 0 | 1.00000 | 0.0136640 | 0.0500307 | 1.0849479 | -3.189938 | 1.9001495 | 5.090087 | -0.4519034 | -0.1114963 | 0.0778499 |
| 4B_ii | 30 | 165 | 0 | 1.00000 | 0.0906950 | 0.0951863 | 1.2085803 | -2.144976 | 1.1713697 | 3.316346 | -0.5656116 | -0.7852813 | 0.0778499 |
| 4B | 31 | 165 | 0 | 1.00000 | 0.0337660 | 0.0538274 | 1.2230776 | -2.091523 | 1.5667310 | 3.658254 | -0.2971274 | -1.0073971 | 0.0778499 |
| 4C | 32 | 165 | 0 | 1.00000 | 0.3156848 | 0.3066383 | 0.0000000 | -3.726638 | 0.3156848 | 4.042323 | -3.0752792 | 7.9229438 | 0.0778499 |
| 4D_i | 33 | 165 | 0 | 1.00000 | 0.0105641 | 0.0206506 | 1.2586268 | -1.937272 | 1.3779341 | 3.315206 | 0.0106530 | -1.3489840 | 0.0778499 |
| 4D_ii | 34 | 164 | 0 | 1.00000 | -0.1593088 | -0.0670358 | 1.2406344 | -1.275037 | 2.3510805 | 3.626118 | 0.4410753 | -0.9999863 | 0.0780869 |
| 4D_iii | 35 | 165 | 0 | 1.00000 | 0.6084056 | 0.1350804 | 0.3684544 | -2.101636 | 0.8569247 | 2.958561 | -0.9762034 | -0.6092913 | 0.0778499 |
| 4D_iv | 36 | 165 | 0 | 1.00000 | -0.0052307 | 0.0312418 | 0.9772973 | -2.728244 | 1.6855707 | 4.413814 | -0.2607805 | -0.5797391 | 0.0778499 |
| 4D | 37 | 165 | 0 | 1.00000 | -0.0020097 | 0.0400835 | 1.1927732 | -2.267040 | 1.7933556 | 4.060396 | -0.2500926 | -0.9104207 | 0.0778499 |
| 5A_i | 38 | 151 | 0 | 1.00000 | 0.1474877 | 0.1623309 | 1.0651116 | -2.726144 | 0.8658956 | 3.592039 | -1.0477351 | 0.1943839 | 0.0813788 |
| 5A_ii | 39 | 162 | 0 | 1.00000 | 0.2117501 | 0.1593764 | 0.8618404 | -2.916245 | 1.0290886 | 3.945334 | -1.3305136 | 1.4050990 | 0.0785674 |
| 5A_iii | 40 | 163 | 0 | 1.00000 | 0.2867713 | 0.1824176 | 0.6537947 | -3.682035 | 0.7277498 | 4.409785 | -1.7118753 | 2.7077792 | 0.0783260 |
| 5A | 41 | 164 | 0 | 1.00000 | 0.1993974 | 0.1282455 | 0.7943910 | -4.746393 | 1.3838408 | 6.130234 | -1.5485784 | 3.4285101 | 0.0780869 |
| 5B_i | 42 | 165 | 0 | 1.00000 | 0.0103107 | 0.0252024 | 1.0297930 | -2.536504 | 2.0199790 | 4.556483 | -0.2416676 | -0.6115278 | 0.0778499 |
| 5B_ii | 43 | 165 | 0 | 1.00000 | 0.0882159 | 0.0241623 | 1.0412917 | -2.865200 | 2.2127576 | 5.077958 | -0.2418312 | -0.2328901 | 0.0778499 |
| 5B_iii | 44 | 141 | 0 | 1.00000 | 0.1195175 | 0.0965819 | 0.8723653 | -3.992784 | 1.7959772 | 5.788761 | -1.1176043 | 2.0818655 | 0.0842152 |
| 5B_iv | 45 | 165 | 0 | 1.00000 | 0.0813121 | 0.1123508 | 1.5300998 | -3.014802 | 1.1133502 | 4.128153 | -0.6830197 | -0.1333095 | 0.0778499 |
| 5B_v | 46 | 162 | 0 | 1.00000 | 0.2069391 | 0.0852912 | 0.9779416 | -2.290919 | 1.2366773 | 3.527596 | -0.6668615 | -0.5273101 | 0.0785674 |
| 5B_vi | 47 | 165 | 0 | 1.00000 | 0.8060507 | 0.0922849 | 0.0000000 | -1.573169 | 0.8060507 | 2.379220 | -0.5393010 | -1.5283369 | 0.0778499 |
| 5B_vii | 48 | 165 | 0 | 1.00000 | 0.0635383 | 0.0637392 | 0.8359201 | -3.356763 | 2.0880643 | 5.444827 | -0.6917272 | 0.6248359 | 0.0778499 |
| 5B | 49 | 165 | 0 | 1.00000 | 0.0745501 | 0.0068070 | 1.1514879 | -2.932793 | 2.3613510 | 5.294145 | -0.1312791 | -0.2710623 | 0.0778499 |
| 5C_i | 50 | 158 | 0 | 1.00000 | 0.0405129 | -0.0245570 | 1.0110083 | -2.461305 | 2.6907121 | 5.152017 | 0.2172977 | -0.0185238 | 0.0795557 |
| 5C_ii | 51 | 165 | 0 | 1.00000 | 0.1343043 | 0.0709744 | 1.0406477 | -2.272235 | 2.0394813 | 4.311716 | -0.5272439 | -0.0966855 | 0.0778499 |
| 5C_iii | 52 | 162 | 0 | 1.00000 | 0.1816585 | 0.0163082 | 1.1694463 | -2.127263 | 1.7300058 | 3.857269 | -0.1625348 | -1.0427012 | 0.0785674 |
| 5C_iv | 53 | 165 | 0 | 1.00000 | 0.1088708 | 0.0363552 | 0.9589692 | -3.067028 | 2.0652579 | 5.132285 | -0.4761717 | 0.7117627 | 0.0778499 |
| 5C | 54 | 165 | 0 | 1.00000 | 0.0364074 | 0.0307456 | 0.9180388 | -2.950263 | 2.5065505 | 5.456813 | -0.2826429 | -0.0730321 | 0.0778499 |
| 5D_i | 55 | 165 | 0 | 1.00000 | 0.0309025 | 0.0401575 | 1.0183139 | -2.688780 | 2.1862212 | 4.875001 | -0.2659083 | -0.3002707 | 0.0778499 |
| 5D_ii | 56 | 161 | 0 | 1.00000 | 0.2099725 | 0.1136750 | 0.8577488 | -2.986832 | 1.4130580 | 4.399890 | -1.0296176 | 0.7555094 | 0.0788110 |
| 5D_iii | 57 | 165 | 0 | 1.00000 | -0.0876865 | 0.0488776 | 0.9149596 | -2.271432 | 2.0960589 | 4.367491 | -0.3781137 | -0.4245462 | 0.0778499 |
| 5D | 58 | 165 | 0 | 1.00000 | 0.0711486 | 0.0531656 | 0.8390709 | -2.731732 | 2.3260600 | 5.057792 | -0.4185826 | 0.0440479 | 0.0778499 |
# export
write_xlsx(x = mydata_std,
path = "economic_freedom_cleaned_standardised.xlsx")# 1. Rename variables to make them lavaan-safe (prepend "X")
mydata_std <- mydata_std %>%
dplyr::rename_with(.fn = ~ paste0("X", .x), .cols = matches("^[0-9]"))
# check: now you have X1A, X1B, X1C, X1D_i, ...
# names(mydata_std)[1:20]Models
CFA
Once the data has been normalized (see previous section), we build the models to test the assumption of equal weighting at three levels:
Across the five Areas (Area 1–Area 5),
Across the subcategories within each Area, and
Across the individual indicators within each subcategory.
Actual (current) Model
- Cannot multiply within
lavaan
mod <- '
Economic_Freedom =~ Area1 + Area2 + Area3 + Area4 + Area5
# Area 1
Area1 =~ X1A + X1B + X1C + X1D + X1E
# Area 2
Area2 =~ (X2A + X2B + X2C + X2D + X2E + X2F + X2G + X2H) * (1 + X2_adj)/2
# Area 3
Area3 =~ X3A + X3B + X3C + X3D
# Area 4
Area4 =~ X4A + X4B + X4C + X4D
# Area 5
Area5 =~ X5A + X5B + X5C + X5D
'
?cfa
fit <- cfa(model = mod,
data = mydata_std,
estimator = "MLR",
std.lv = TRUE,
missing = "fiml"
)
summary(fit, standardized = TRUE, fit.measures = TRUE)
# --- (Optional) factor scores for Areas and EFI ---
scores <- lavPredict(fit) # matrix with Area1..Area5 and EFI scores
head(scores)
# --- (Optional) path diagram ---
lavaanPlot(model = fit, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))So I cannot run this model. Lets try to approximate it.
Closest Possible Model?
Keeping x2_adj as a variable for
Area 2.
mod <- '
Economic_Freedom =~ Area1 + Area2 + Area3 + Area4 + Area5
# Area 1
Area1 =~ X1A + X1B + X1C + X1D + X1E
# Area 2
Area2 =~ X2A + X2B + X2C + X2D + X2E + X2F + X2G + X2H + X2_adj
# Area 3
Area3 =~ X3A + X3B + X3C + X3D
# Area 4
Area4 =~ X4A + X4B + X4C + X4D
# Area 5
Area5 =~ X5A + X5B + X5C + X5D
'
?cfa## Help on topic 'cfa' was found in the following packages:
##
## Package Library
## sem /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## lavaan /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
##
##
## Using the first match ...fit <- cfa(model = mod,
data = mydata_std,
estimator = "MLR",
std.lv = TRUE,
missing = "fiml"
)
summary(fit, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 110 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 83
##
## Number of observations 165
## Number of missing patterns 13
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 1133.244 1027.213
## Degrees of freedom 294 294
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.103
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 3473.006 3108.849
## Degrees of freedom 325 325
## P-value 0.000 0.000
## Scaling correction factor 1.117
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.733 0.737
## Tucker-Lewis Index (TLI) 0.705 0.709
##
## Robust Comparative Fit Index (CFI) 0.740
## Robust Tucker-Lewis Index (TLI) 0.712
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -4834.797 -4834.797
## Scaling correction factor 1.208
## for the MLR correction
## Loglikelihood unrestricted model (H1) -4268.175 -4268.175
## Scaling correction factor 1.126
## for the MLR correction
##
## Akaike (AIC) 9835.594 9835.594
## Bayesian (BIC) 10093.388 10093.388
## Sample-size adjusted Bayesian (SABIC) 9830.609 9830.609
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.132 0.123
## 90 Percent confidence interval - lower 0.123 0.115
## 90 Percent confidence interval - upper 0.140 0.131
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.130
## 90 Percent confidence interval - lower 0.121
## 90 Percent confidence interval - upper 0.139
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.118 0.118
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 -1.648 0.607 -2.716 0.007 -0.855 -0.855
## Area2 6.230 3.885 1.604 0.109 0.987 0.987
## Area3 0.429 0.113 3.807 0.000 0.395 0.395
## Area4 1.957 0.370 5.294 0.000 0.890 0.890
## Area5 4.073 1.163 3.503 0.000 0.971 0.971
## Area1 =~
## X1A 0.269 0.072 3.738 0.000 0.519 0.520
## X1B 0.369 0.141 2.620 0.009 0.710 0.715
## X1C -0.232 0.075 -3.085 0.002 -0.447 -0.449
## X1D 0.184 0.091 2.019 0.043 0.355 0.356
## X1E -0.218 0.067 -3.253 0.001 -0.420 -0.422
## Area2 =~
## X2A 0.142 0.088 1.616 0.106 0.897 0.900
## X2B 0.151 0.092 1.640 0.101 0.953 0.956
## X2C 0.144 0.087 1.655 0.098 0.907 0.910
## X2D 0.122 0.073 1.669 0.095 0.771 0.788
## X2E 0.142 0.085 1.670 0.095 0.894 0.899
## X2F 0.130 0.078 1.660 0.097 0.820 0.822
## X2G 0.050 0.032 1.577 0.115 0.314 0.315
## X2H 0.124 0.075 1.659 0.097 0.786 0.788
## X2_adj 0.071 0.043 1.657 0.097 0.447 0.448
## Area3 =~
## X3A 0.650 0.156 4.170 0.000 0.707 0.711
## X3B 0.462 0.113 4.088 0.000 0.502 0.504
## X3C 0.755 0.071 10.612 0.000 0.822 0.825
## X3D 0.104 0.082 1.273 0.203 0.113 0.113
## Area4 =~
## X4A 0.223 0.051 4.408 0.000 0.490 0.492
## X4B 0.399 0.066 6.008 0.000 0.877 0.880
## X4C 0.172 0.049 3.523 0.000 0.378 0.379
## X4D 0.372 0.052 7.123 0.000 0.817 0.820
## Area5 =~
## X5A 0.098 0.033 2.980 0.003 0.410 0.411
## X5B 0.124 0.035 3.550 0.000 0.521 0.522
## X5C 0.221 0.062 3.579 0.000 0.927 0.930
## X5D 0.189 0.051 3.687 0.000 0.791 0.793
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1B 0.019 0.078 0.242 0.809 0.019 0.019
## .X1C -0.005 0.079 -0.065 0.948 -0.005 -0.005
## .X1D -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1E 0.000 0.078 0.002 0.998 0.000 0.000
## .X2A 0.000 0.078 0.000 1.000 0.000 0.000
## .X2B 0.000 0.078 0.000 1.000 0.000 0.000
## .X2C 0.000 0.078 0.000 1.000 0.000 0.000
## .X2D -0.059 0.080 -0.745 0.456 -0.059 -0.061
## .X2E -0.002 0.078 -0.032 0.974 -0.002 -0.002
## .X2F 0.000 0.078 0.000 1.000 0.000 0.000
## .X2G -0.003 0.078 -0.043 0.966 -0.003 -0.003
## .X2H -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2_adj 0.000 0.078 0.000 1.000 0.000 0.000
## .X3A 0.004 0.077 0.048 0.962 0.004 0.004
## .X3B 0.000 0.078 0.000 1.000 0.000 0.000
## .X3C 0.000 0.078 0.000 1.000 0.000 0.000
## .X3D 0.000 0.078 0.000 1.000 0.000 0.000
## .X4A 0.000 0.078 0.000 1.000 0.000 0.000
## .X4B 0.000 0.078 0.000 1.000 0.000 0.000
## .X4C 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D 0.000 0.078 0.000 1.000 0.000 0.000
## .X5A -0.005 0.078 -0.066 0.947 -0.005 -0.005
## .X5B 0.000 0.078 0.000 1.000 0.000 0.000
## .X5C 0.000 0.078 0.000 1.000 0.000 0.000
## .X5D -0.000 0.078 -0.000 1.000 -0.000 -0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.725 0.114 6.385 0.000 0.725 0.729
## .X1B 0.483 0.122 3.944 0.000 0.483 0.489
## .X1C 0.793 0.093 8.519 0.000 0.793 0.799
## .X1D 0.868 0.085 10.159 0.000 0.868 0.873
## .X1E 0.816 0.096 8.492 0.000 0.816 0.822
## .X2A 0.189 0.026 7.390 0.000 0.189 0.190
## .X2B 0.086 0.018 4.878 0.000 0.086 0.087
## .X2C 0.171 0.031 5.420 0.000 0.171 0.172
## .X2D 0.363 0.042 8.594 0.000 0.363 0.379
## .X2E 0.191 0.026 7.366 0.000 0.191 0.192
## .X2F 0.322 0.042 7.760 0.000 0.322 0.324
## .X2G 0.896 0.106 8.424 0.000 0.896 0.901
## .X2H 0.377 0.051 7.338 0.000 0.377 0.379
## .X2_adj 0.794 0.116 6.847 0.000 0.794 0.799
## .X3A 0.491 0.139 3.518 0.000 0.491 0.495
## .X3B 0.742 0.149 4.966 0.000 0.742 0.746
## .X3C 0.318 0.103 3.102 0.002 0.318 0.320
## .X3D 0.981 0.091 10.803 0.000 0.981 0.987
## .X4A 0.754 0.102 7.408 0.000 0.754 0.758
## .X4B 0.225 0.051 4.426 0.000 0.225 0.226
## .X4C 0.851 0.183 4.655 0.000 0.851 0.857
## .X4D 0.326 0.055 5.933 0.000 0.326 0.328
## .X5A 0.829 0.132 6.259 0.000 0.829 0.831
## .X5B 0.723 0.078 9.256 0.000 0.723 0.727
## .X5C 0.134 0.025 5.441 0.000 0.134 0.135
## .X5D 0.368 0.056 6.539 0.000 0.368 0.371
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.269 0.269
## .Area2 1.000 0.025 0.025
## .Area3 1.000 0.844 0.844
## .Area4 1.000 0.207 0.207
## .Area5 1.000 0.057 0.057# --- (Optional) factor scores for Areas and EFI ---
scores <- lavPredict(fit) # matrix with Area1..Area5 and EFI scores
head(scores)## Economic_Freedom Area1 Area2 Area3 Area4 Area5
## [1,] 0.05901158 0.0203957 -0.1427195 0.64637879 1.4440979 0.3792955
## [2,] -0.94707730 1.3954745 -5.3714503 -0.06755297 -3.7185180 -4.1049904
## [3,] -1.20192036 2.3689464 -7.4228897 -1.41822796 -3.4916190 -4.4897822
## [4,] -0.40937561 -0.2319512 -2.2480537 -3.28109534 -1.6386597 -1.8672775
## [5,] 0.08875661 -0.0369960 0.1902693 0.31269536 0.7601009 0.6715767
## [6,] 1.73738792 -2.8202697 11.4863358 0.73512001 2.4271128 6.9747549# --- (Optional) path diagram ---
lavaanPlot(model = fit, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))Lowest Level Model
- All individual indicators within each subcategory.
mod <- '
Economic_Freedom =~ Area1 + Area2 + Area3 + Area4 + Area5
# Area 1
Area1 =~ X1A + X1B + X1C + X1D_i + X1D_ii + X1E
# Area 2
Area2 =~ X2A + X2B + X2C + X2D + X2E + X2F + X2G + X2H + X2_adj
# Area 3
Area3 =~ X3A + X3B + X3C + X3D
# Area 4
Area4 =~ X4A_i + X4A_ii + X4A_iii +
X4B_i + X4B_ii +
X4C +
X4D_i + X4D_ii + X4D_iii + X4D_iv
# Area 5
Area5 =~ X5A_i + X5A_ii + X5A_iii +
X5B_i + X5B_ii + X5B_iii + X5B_iv + X5B_v + X5B_vi + X5B_vii +
X5C_i + X5C_ii + X5C_iii + X5C_iv +
X5D_i + X5D_ii + X5D_iii
'
fit <- cfa(mod, data = mydata_std, estimator = "MLR", std.lv = TRUE, missing = "fiml")
summary(fit, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 126 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 143
##
## Number of observations 165
## Number of missing patterns 31
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 3166.420 3008.129
## Degrees of freedom 984 984
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.053
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 6413.070 6013.842
## Degrees of freedom 1035 1035
## P-value 0.000 0.000
## Scaling correction factor 1.066
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.594 0.593
## Tucker-Lewis Index (TLI) 0.573 0.572
##
## Robust Comparative Fit Index (CFI) 0.599
## Robust Tucker-Lewis Index (TLI) 0.579
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -8950.237 -8950.237
## Scaling correction factor 1.200
## for the MLR correction
## Loglikelihood unrestricted model (H1) -7367.027 -7367.027
## Scaling correction factor 1.071
## for the MLR correction
##
## Akaike (AIC) 18186.474 18186.474
## Bayesian (BIC) 18630.624 18630.624
## Sample-size adjusted Bayesian (SABIC) 18177.885 18177.885
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.116 0.112
## 90 Percent confidence interval - lower 0.111 0.107
## 90 Percent confidence interval - upper 0.120 0.116
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.115
## 90 Percent confidence interval - lower 0.110
## 90 Percent confidence interval - upper 0.120
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.114 0.114
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 -1.195 0.342 -3.493 0.000 -0.767 -0.767
## Area2 2.980 0.625 4.769 0.000 0.948 0.948
## Area3 0.449 0.120 3.741 0.000 0.410 0.410
## Area4 2.365 0.492 4.805 0.000 0.921 0.921
## Area5 4.909 2.037 2.410 0.016 0.980 0.980
## Area1 =~
## X1A 0.287 0.062 4.634 0.000 0.448 0.449
## X1B 0.485 0.109 4.463 0.000 0.756 0.761
## X1C -0.295 0.063 -4.651 0.000 -0.459 -0.461
## X1D_i 0.301 0.113 2.666 0.008 0.469 0.470
## X1D_ii 0.273 0.120 2.270 0.023 0.426 0.427
## X1E -0.288 0.066 -4.351 0.000 -0.449 -0.450
## Area2 =~
## X2A 0.286 0.054 5.344 0.000 0.899 0.902
## X2B 0.304 0.055 5.562 0.000 0.956 0.959
## X2C 0.290 0.050 5.783 0.000 0.911 0.914
## X2D 0.244 0.044 5.595 0.000 0.768 0.783
## X2E 0.283 0.050 5.653 0.000 0.889 0.893
## X2F 0.261 0.047 5.600 0.000 0.820 0.823
## X2G 0.098 0.028 3.475 0.001 0.307 0.308
## X2H 0.248 0.044 5.688 0.000 0.780 0.782
## X2_adj 0.141 0.030 4.757 0.000 0.442 0.444
## Area3 =~
## X3A 0.647 0.153 4.217 0.000 0.709 0.712
## X3B 0.463 0.111 4.186 0.000 0.508 0.509
## X3C 0.744 0.069 10.752 0.000 0.815 0.818
## X3D 0.109 0.081 1.345 0.179 0.120 0.120
## Area4 =~
## X4A_i 0.220 0.052 4.276 0.000 0.566 0.567
## X4A_ii 0.256 0.056 4.536 0.000 0.657 0.659
## X4A_iii 0.063 0.034 1.830 0.067 0.162 0.163
## X4B_i 0.317 0.061 5.196 0.000 0.813 0.816
## X4B_ii 0.269 0.046 5.915 0.000 0.691 0.694
## X4C 0.143 0.042 3.432 0.001 0.368 0.369
## X4D_i 0.346 0.064 5.424 0.000 0.888 0.891
## X4D_ii 0.197 0.050 3.968 0.000 0.506 0.507
## X4D_iii 0.152 0.038 4.027 0.000 0.391 0.392
## X4D_iv 0.341 0.056 6.114 0.000 0.876 0.879
## Area5 =~
## X5A_i 0.073 0.033 2.204 0.028 0.363 0.364
## X5A_ii 0.040 0.021 1.899 0.058 0.198 0.199
## X5A_iii 0.037 0.023 1.564 0.118 0.183 0.184
## X5B_i 0.080 0.041 1.944 0.052 0.401 0.402
## X5B_ii 0.104 0.050 2.060 0.039 0.519 0.521
## X5B_iii 0.026 0.025 1.030 0.303 0.128 0.128
## X5B_iv 0.019 0.019 1.002 0.316 0.094 0.095
## X5B_v 0.070 0.034 2.054 0.040 0.350 0.350
## X5B_vi 0.026 0.019 1.362 0.173 0.132 0.133
## X5B_vii 0.122 0.053 2.311 0.021 0.611 0.613
## X5C_i 0.079 0.036 2.221 0.026 0.397 0.398
## X5C_ii 0.162 0.063 2.588 0.010 0.810 0.813
## X5C_iii 0.148 0.055 2.715 0.007 0.743 0.748
## X5C_iv 0.137 0.058 2.371 0.018 0.688 0.690
## X5D_i 0.168 0.071 2.375 0.018 0.843 0.845
## X5D_ii 0.041 0.028 1.490 0.136 0.205 0.206
## X5D_iii 0.151 0.066 2.277 0.023 0.758 0.760
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1B 0.018 0.078 0.235 0.814 0.018 0.019
## .X1C -0.006 0.079 -0.074 0.941 -0.006 -0.006
## .X1D_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1D_ii 0.007 0.079 0.086 0.932 0.007 0.007
## .X1E -0.002 0.079 -0.024 0.981 -0.002 -0.002
## .X2A 0.000 0.078 0.001 1.000 0.000 0.000
## .X2B 0.000 0.078 0.001 0.999 0.000 0.000
## .X2C 0.000 0.078 0.001 1.000 0.000 0.000
## .X2D -0.062 0.080 -0.782 0.434 -0.062 -0.064
## .X2E -0.003 0.078 -0.035 0.972 -0.003 -0.003
## .X2F 0.000 0.078 0.001 1.000 0.000 0.000
## .X2G -0.003 0.078 -0.042 0.966 -0.003 -0.003
## .X2H 0.000 0.078 0.001 1.000 0.000 0.000
## .X2_adj 0.000 0.078 0.000 1.000 0.000 0.000
## .X3A 0.004 0.077 0.049 0.961 0.004 0.004
## .X3B 0.000 0.078 0.000 1.000 0.000 0.000
## .X3C 0.000 0.078 0.000 1.000 0.000 0.000
## .X3D 0.000 0.078 0.000 1.000 0.000 0.000
## .X4A_i -0.024 0.081 -0.299 0.765 -0.024 -0.024
## .X4A_ii -0.001 0.078 -0.018 0.986 -0.001 -0.001
## .X4A_iii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4B_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X4B_ii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4C 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X4D_ii 0.003 0.078 0.039 0.969 0.003 0.003
## .X4D_iii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D_iv 0.000 0.078 0.001 1.000 0.000 0.000
## .X5A_i -0.029 0.084 -0.342 0.732 -0.029 -0.029
## .X5A_ii -0.005 0.079 -0.061 0.952 -0.005 -0.005
## .X5A_iii -0.004 0.079 -0.048 0.962 -0.004 -0.004
## .X5B_i 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_ii 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_iii -0.017 0.080 -0.211 0.833 -0.017 -0.017
## .X5B_iv 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_v -0.007 0.079 -0.094 0.925 -0.007 -0.007
## .X5B_vi 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_vii 0.000 0.078 0.000 1.000 0.000 0.000
## .X5C_i -0.019 0.079 -0.246 0.806 -0.019 -0.019
## .X5C_ii 0.000 0.078 0.001 1.000 0.000 0.000
## .X5C_iii 0.003 0.078 0.037 0.970 0.003 0.003
## .X5C_iv 0.000 0.078 0.000 1.000 0.000 0.000
## .X5D_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X5D_ii -0.008 0.079 -0.107 0.915 -0.008 -0.008
## .X5D_iii 0.000 0.078 0.001 1.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.793 0.126 6.304 0.000 0.793 0.798
## .X1B 0.415 0.099 4.197 0.000 0.415 0.421
## .X1C 0.781 0.089 8.813 0.000 0.781 0.787
## .X1D_i 0.774 0.101 7.670 0.000 0.774 0.779
## .X1D_ii 0.813 0.096 8.429 0.000 0.813 0.818
## .X1E 0.791 0.094 8.377 0.000 0.791 0.797
## .X2A 0.186 0.025 7.313 0.000 0.186 0.187
## .X2B 0.081 0.018 4.569 0.000 0.081 0.081
## .X2C 0.164 0.033 4.991 0.000 0.164 0.165
## .X2D 0.371 0.043 8.583 0.000 0.371 0.386
## .X2E 0.200 0.027 7.428 0.000 0.200 0.202
## .X2F 0.321 0.042 7.651 0.000 0.321 0.323
## .X2G 0.899 0.107 8.404 0.000 0.899 0.905
## .X2H 0.386 0.053 7.257 0.000 0.386 0.388
## .X2_adj 0.798 0.116 6.863 0.000 0.798 0.803
## .X3A 0.488 0.138 3.528 0.000 0.488 0.492
## .X3B 0.736 0.149 4.935 0.000 0.736 0.740
## .X3C 0.329 0.096 3.413 0.001 0.329 0.331
## .X3D 0.980 0.091 10.736 0.000 0.980 0.986
## .X4A_i 0.675 0.121 5.578 0.000 0.675 0.678
## .X4A_ii 0.561 0.116 4.814 0.000 0.561 0.565
## .X4A_iii 0.968 0.154 6.301 0.000 0.968 0.974
## .X4B_i 0.333 0.058 5.770 0.000 0.333 0.335
## .X4B_ii 0.516 0.068 7.586 0.000 0.516 0.519
## .X4C 0.858 0.186 4.610 0.000 0.858 0.864
## .X4D_i 0.206 0.041 5.075 0.000 0.206 0.207
## .X4D_ii 0.738 0.103 7.181 0.000 0.738 0.743
## .X4D_iii 0.841 0.084 9.985 0.000 0.841 0.846
## .X4D_iv 0.226 0.040 5.585 0.000 0.226 0.227
## .X5A_i 0.866 0.093 9.318 0.000 0.866 0.868
## .X5A_ii 0.955 0.138 6.923 0.000 0.955 0.961
## .X5A_iii 0.961 0.152 6.323 0.000 0.961 0.966
## .X5B_i 0.833 0.086 9.672 0.000 0.833 0.838
## .X5B_ii 0.725 0.082 8.875 0.000 0.725 0.729
## .X5B_iii 0.977 0.178 5.481 0.000 0.977 0.984
## .X5B_iv 0.985 0.106 9.270 0.000 0.985 0.991
## .X5B_v 0.876 0.089 9.803 0.000 0.876 0.878
## .X5B_vi 0.976 0.056 17.370 0.000 0.976 0.982
## .X5B_vii 0.621 0.086 7.244 0.000 0.621 0.625
## .X5C_i 0.841 0.093 9.045 0.000 0.841 0.842
## .X5C_ii 0.337 0.039 8.604 0.000 0.337 0.339
## .X5C_iii 0.434 0.072 6.051 0.000 0.434 0.440
## .X5C_iv 0.521 0.068 7.649 0.000 0.521 0.524
## .X5D_i 0.284 0.042 6.681 0.000 0.284 0.285
## .X5D_ii 0.954 0.127 7.503 0.000 0.954 0.958
## .X5D_iii 0.420 0.071 5.913 0.000 0.420 0.422
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.412 0.412
## .Area2 1.000 0.101 0.101
## .Area3 1.000 0.832 0.832
## .Area4 1.000 0.152 0.152
## .Area5 1.000 0.040 0.040# --- (Optional) factor scores for Areas and EFI ---
scores <- lavPredict(fit) # matrix with Area1..Area5 and EFI scores
head(scores)## Economic_Freedom Area1 Area2 Area3 Area4 Area5
## [1,] 0.2487350 -0.02086186 -0.13771047 0.6728710 1.1842258 1.5339617
## [2,] -1.0796441 0.92747919 -2.54620172 -0.0968184 -3.7245216 -5.4865630
## [3,] -1.2293312 1.91909606 -3.63809460 -1.4199126 -3.3778351 -5.8846831
## [4,] -0.5462043 -0.59492358 -1.08571834 -3.3303559 -2.5060237 -2.5579467
## [5,] 0.1057783 0.09580203 -0.03691889 0.3185355 0.4824951 0.6719004
## [6,] 1.5933738 -2.00952653 5.81610096 0.7336082 3.0544472 7.8146341# --- (Optional) path diagram ---
lavaanPlot(model = fit, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))Free vs Equal Model
Ideally, we would like to check/show that the equal model of EFI is significantly worse than free model.
Attempt 1
# -------- FREE ----------
mod_free <- '
Economic_Freedom =~ Area1 + Area2 + Area3 + Area4 + Area5
Area1 =~ X1A + X1B + X1C + X1D_i + X1D_ii + X1E
Area2 =~ X2A + X2B + X2C + X2D + X2E + X2F + X2G + X2H + X2_adj
Area3 =~ X3A + X3B + X3C + X3D
Area4 =~ X4A_i + X4A_ii + X4A_iii + X4B_i + X4B_ii + X4C +
X4D_i + X4D_ii + X4D_iii + X4D_iv
Area5 =~ X5A_i + X5A_ii + X5A_iii +
X5B_i + X5B_ii + X5B_iii + X5B_iv + X5B_v + X5B_vi + X5B_vii +
X5C_i + X5C_ii + X5C_iii + X5C_iv +
X5D_i + X5D_ii + X5D_iii
'
# -------- EQUAL ----------
mod_equal <- '
# Top level: each Area contributes equally to Economic_Freedom
Economic_Freedom =~ 1*Area1 + 1*Area2 + 1*Area3 + 1*Area4 + 1*Area5
# --- Area 1: 5 indicators → each weight = 1/5 ---
Area1 =~ 0.2*X1A + 0.2*X1B + 0.2*X1C + 0.2*X1D + 0.2*X1E
# --- Area 2: 8 indicators → each weight = 1/8 ---
Area2 =~ 0.125*X2A + 0.125*X2B + 0.125*X2C + 0.125*X2D +
0.125*X2E + 0.125*X2F + 0.125*X2G + 0.125*X2H + X2_adj
# --- Area 3: 4 indicators → each weight = 1/4 ---
Area3 =~ 0.25*X3A + 0.25*X3B + 0.25*X3C + 0.25*X3D
# --- Area 4: 4 indicators → each weight = 1/4 ---
Area4 =~ 0.25*X4A + 0.25*X4B + 0.25*X4C + 0.25*X4D
# --- Area 5: 4 indicators → each weight = 1/4 ---
Area5 =~ 0.25*X5A + 0.25*X5B + 0.25*X5C + 0.25*X5D
'
# -------- Fit & compare ----------
dat <- mydata_std
fit_free <- cfa(mod_free, data = dat, estimator = "MLR", std.lv = TRUE, missing = "fiml")
fit_equal <- cfa(mod_equal, data = dat, estimator = "MLR", std.lv = TRUE, missing = "fiml")
summary(fit_free, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 126 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 143
##
## Number of observations 165
## Number of missing patterns 31
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 3166.420 3008.129
## Degrees of freedom 984 984
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.053
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 6413.070 6013.842
## Degrees of freedom 1035 1035
## P-value 0.000 0.000
## Scaling correction factor 1.066
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.594 0.593
## Tucker-Lewis Index (TLI) 0.573 0.572
##
## Robust Comparative Fit Index (CFI) 0.599
## Robust Tucker-Lewis Index (TLI) 0.579
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -8950.237 -8950.237
## Scaling correction factor 1.200
## for the MLR correction
## Loglikelihood unrestricted model (H1) -7367.027 -7367.027
## Scaling correction factor 1.071
## for the MLR correction
##
## Akaike (AIC) 18186.474 18186.474
## Bayesian (BIC) 18630.624 18630.624
## Sample-size adjusted Bayesian (SABIC) 18177.885 18177.885
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.116 0.112
## 90 Percent confidence interval - lower 0.111 0.107
## 90 Percent confidence interval - upper 0.120 0.116
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.115
## 90 Percent confidence interval - lower 0.110
## 90 Percent confidence interval - upper 0.120
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.114 0.114
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 -1.195 0.342 -3.493 0.000 -0.767 -0.767
## Area2 2.980 0.625 4.769 0.000 0.948 0.948
## Area3 0.449 0.120 3.741 0.000 0.410 0.410
## Area4 2.365 0.492 4.805 0.000 0.921 0.921
## Area5 4.909 2.037 2.410 0.016 0.980 0.980
## Area1 =~
## X1A 0.287 0.062 4.634 0.000 0.448 0.449
## X1B 0.485 0.109 4.463 0.000 0.756 0.761
## X1C -0.295 0.063 -4.651 0.000 -0.459 -0.461
## X1D_i 0.301 0.113 2.666 0.008 0.469 0.470
## X1D_ii 0.273 0.120 2.270 0.023 0.426 0.427
## X1E -0.288 0.066 -4.351 0.000 -0.449 -0.450
## Area2 =~
## X2A 0.286 0.054 5.344 0.000 0.899 0.902
## X2B 0.304 0.055 5.562 0.000 0.956 0.959
## X2C 0.290 0.050 5.783 0.000 0.911 0.914
## X2D 0.244 0.044 5.595 0.000 0.768 0.783
## X2E 0.283 0.050 5.653 0.000 0.889 0.893
## X2F 0.261 0.047 5.600 0.000 0.820 0.823
## X2G 0.098 0.028 3.475 0.001 0.307 0.308
## X2H 0.248 0.044 5.688 0.000 0.780 0.782
## X2_adj 0.141 0.030 4.757 0.000 0.442 0.444
## Area3 =~
## X3A 0.647 0.153 4.217 0.000 0.709 0.712
## X3B 0.463 0.111 4.186 0.000 0.508 0.509
## X3C 0.744 0.069 10.752 0.000 0.815 0.818
## X3D 0.109 0.081 1.345 0.179 0.120 0.120
## Area4 =~
## X4A_i 0.220 0.052 4.276 0.000 0.566 0.567
## X4A_ii 0.256 0.056 4.536 0.000 0.657 0.659
## X4A_iii 0.063 0.034 1.830 0.067 0.162 0.163
## X4B_i 0.317 0.061 5.196 0.000 0.813 0.816
## X4B_ii 0.269 0.046 5.915 0.000 0.691 0.694
## X4C 0.143 0.042 3.432 0.001 0.368 0.369
## X4D_i 0.346 0.064 5.424 0.000 0.888 0.891
## X4D_ii 0.197 0.050 3.968 0.000 0.506 0.507
## X4D_iii 0.152 0.038 4.027 0.000 0.391 0.392
## X4D_iv 0.341 0.056 6.114 0.000 0.876 0.879
## Area5 =~
## X5A_i 0.073 0.033 2.204 0.028 0.363 0.364
## X5A_ii 0.040 0.021 1.899 0.058 0.198 0.199
## X5A_iii 0.037 0.023 1.564 0.118 0.183 0.184
## X5B_i 0.080 0.041 1.944 0.052 0.401 0.402
## X5B_ii 0.104 0.050 2.060 0.039 0.519 0.521
## X5B_iii 0.026 0.025 1.030 0.303 0.128 0.128
## X5B_iv 0.019 0.019 1.002 0.316 0.094 0.095
## X5B_v 0.070 0.034 2.054 0.040 0.350 0.350
## X5B_vi 0.026 0.019 1.362 0.173 0.132 0.133
## X5B_vii 0.122 0.053 2.311 0.021 0.611 0.613
## X5C_i 0.079 0.036 2.221 0.026 0.397 0.398
## X5C_ii 0.162 0.063 2.588 0.010 0.810 0.813
## X5C_iii 0.148 0.055 2.715 0.007 0.743 0.748
## X5C_iv 0.137 0.058 2.371 0.018 0.688 0.690
## X5D_i 0.168 0.071 2.375 0.018 0.843 0.845
## X5D_ii 0.041 0.028 1.490 0.136 0.205 0.206
## X5D_iii 0.151 0.066 2.277 0.023 0.758 0.760
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1B 0.018 0.078 0.235 0.814 0.018 0.019
## .X1C -0.006 0.079 -0.074 0.941 -0.006 -0.006
## .X1D_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X1D_ii 0.007 0.079 0.086 0.932 0.007 0.007
## .X1E -0.002 0.079 -0.024 0.981 -0.002 -0.002
## .X2A 0.000 0.078 0.001 1.000 0.000 0.000
## .X2B 0.000 0.078 0.001 0.999 0.000 0.000
## .X2C 0.000 0.078 0.001 1.000 0.000 0.000
## .X2D -0.062 0.080 -0.782 0.434 -0.062 -0.064
## .X2E -0.003 0.078 -0.035 0.972 -0.003 -0.003
## .X2F 0.000 0.078 0.001 1.000 0.000 0.000
## .X2G -0.003 0.078 -0.042 0.966 -0.003 -0.003
## .X2H 0.000 0.078 0.001 1.000 0.000 0.000
## .X2_adj 0.000 0.078 0.000 1.000 0.000 0.000
## .X3A 0.004 0.077 0.049 0.961 0.004 0.004
## .X3B 0.000 0.078 0.000 1.000 0.000 0.000
## .X3C 0.000 0.078 0.000 1.000 0.000 0.000
## .X3D 0.000 0.078 0.000 1.000 0.000 0.000
## .X4A_i -0.024 0.081 -0.299 0.765 -0.024 -0.024
## .X4A_ii -0.001 0.078 -0.018 0.986 -0.001 -0.001
## .X4A_iii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4B_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X4B_ii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4C 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X4D_ii 0.003 0.078 0.039 0.969 0.003 0.003
## .X4D_iii 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D_iv 0.000 0.078 0.001 1.000 0.000 0.000
## .X5A_i -0.029 0.084 -0.342 0.732 -0.029 -0.029
## .X5A_ii -0.005 0.079 -0.061 0.952 -0.005 -0.005
## .X5A_iii -0.004 0.079 -0.048 0.962 -0.004 -0.004
## .X5B_i 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_ii 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_iii -0.017 0.080 -0.211 0.833 -0.017 -0.017
## .X5B_iv 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_v -0.007 0.079 -0.094 0.925 -0.007 -0.007
## .X5B_vi 0.000 0.078 0.000 1.000 0.000 0.000
## .X5B_vii 0.000 0.078 0.000 1.000 0.000 0.000
## .X5C_i -0.019 0.079 -0.246 0.806 -0.019 -0.019
## .X5C_ii 0.000 0.078 0.001 1.000 0.000 0.000
## .X5C_iii 0.003 0.078 0.037 0.970 0.003 0.003
## .X5C_iv 0.000 0.078 0.000 1.000 0.000 0.000
## .X5D_i 0.000 0.078 0.001 1.000 0.000 0.000
## .X5D_ii -0.008 0.079 -0.107 0.915 -0.008 -0.008
## .X5D_iii 0.000 0.078 0.001 1.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.793 0.126 6.304 0.000 0.793 0.798
## .X1B 0.415 0.099 4.197 0.000 0.415 0.421
## .X1C 0.781 0.089 8.813 0.000 0.781 0.787
## .X1D_i 0.774 0.101 7.670 0.000 0.774 0.779
## .X1D_ii 0.813 0.096 8.429 0.000 0.813 0.818
## .X1E 0.791 0.094 8.377 0.000 0.791 0.797
## .X2A 0.186 0.025 7.313 0.000 0.186 0.187
## .X2B 0.081 0.018 4.569 0.000 0.081 0.081
## .X2C 0.164 0.033 4.991 0.000 0.164 0.165
## .X2D 0.371 0.043 8.583 0.000 0.371 0.386
## .X2E 0.200 0.027 7.428 0.000 0.200 0.202
## .X2F 0.321 0.042 7.651 0.000 0.321 0.323
## .X2G 0.899 0.107 8.404 0.000 0.899 0.905
## .X2H 0.386 0.053 7.257 0.000 0.386 0.388
## .X2_adj 0.798 0.116 6.863 0.000 0.798 0.803
## .X3A 0.488 0.138 3.528 0.000 0.488 0.492
## .X3B 0.736 0.149 4.935 0.000 0.736 0.740
## .X3C 0.329 0.096 3.413 0.001 0.329 0.331
## .X3D 0.980 0.091 10.736 0.000 0.980 0.986
## .X4A_i 0.675 0.121 5.578 0.000 0.675 0.678
## .X4A_ii 0.561 0.116 4.814 0.000 0.561 0.565
## .X4A_iii 0.968 0.154 6.301 0.000 0.968 0.974
## .X4B_i 0.333 0.058 5.770 0.000 0.333 0.335
## .X4B_ii 0.516 0.068 7.586 0.000 0.516 0.519
## .X4C 0.858 0.186 4.610 0.000 0.858 0.864
## .X4D_i 0.206 0.041 5.075 0.000 0.206 0.207
## .X4D_ii 0.738 0.103 7.181 0.000 0.738 0.743
## .X4D_iii 0.841 0.084 9.985 0.000 0.841 0.846
## .X4D_iv 0.226 0.040 5.585 0.000 0.226 0.227
## .X5A_i 0.866 0.093 9.318 0.000 0.866 0.868
## .X5A_ii 0.955 0.138 6.923 0.000 0.955 0.961
## .X5A_iii 0.961 0.152 6.323 0.000 0.961 0.966
## .X5B_i 0.833 0.086 9.672 0.000 0.833 0.838
## .X5B_ii 0.725 0.082 8.875 0.000 0.725 0.729
## .X5B_iii 0.977 0.178 5.481 0.000 0.977 0.984
## .X5B_iv 0.985 0.106 9.270 0.000 0.985 0.991
## .X5B_v 0.876 0.089 9.803 0.000 0.876 0.878
## .X5B_vi 0.976 0.056 17.370 0.000 0.976 0.982
## .X5B_vii 0.621 0.086 7.244 0.000 0.621 0.625
## .X5C_i 0.841 0.093 9.045 0.000 0.841 0.842
## .X5C_ii 0.337 0.039 8.604 0.000 0.337 0.339
## .X5C_iii 0.434 0.072 6.051 0.000 0.434 0.440
## .X5C_iv 0.521 0.068 7.649 0.000 0.521 0.524
## .X5D_i 0.284 0.042 6.681 0.000 0.284 0.285
## .X5D_ii 0.954 0.127 7.503 0.000 0.954 0.958
## .X5D_iii 0.420 0.071 5.913 0.000 0.420 0.422
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.412 0.412
## .Area2 1.000 0.101 0.101
## .Area3 1.000 0.832 0.832
## .Area4 1.000 0.152 0.152
## .Area5 1.000 0.040 0.040summary(fit_equal, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 18 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 53
##
## Number of observations 165
## Number of missing patterns 13
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2712.420 2428.214
## Degrees of freedom 324 324
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.117
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 3473.006 3108.849
## Degrees of freedom 325 325
## P-value 0.000 0.000
## Scaling correction factor 1.117
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.241 0.244
## Tucker-Lewis Index (TLI) 0.239 0.242
##
## Robust Comparative Fit Index (CFI) 0.243
## Robust Tucker-Lewis Index (TLI) 0.240
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5624.385 -5624.385
## Scaling correction factor 1.183
## for the MLR correction
## Loglikelihood unrestricted model (H1) -4268.175 -4268.175
## Scaling correction factor 1.126
## for the MLR correction
##
## Akaike (AIC) 11354.770 11354.770
## Bayesian (BIC) 11519.386 11519.386
## Sample-size adjusted Bayesian (SABIC) 11351.587 11351.587
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.211 0.198
## 90 Percent confidence interval - lower 0.204 0.191
## 90 Percent confidence interval - upper 0.219 0.205
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.211
## 90 Percent confidence interval - lower 0.203
## 90 Percent confidence interval - upper 0.219
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.365 0.365
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 1.000 0.707 0.707
## Area2 1.000 0.707 0.707
## Area3 1.000 0.707 0.707
## Area4 1.000 0.707 0.707
## Area5 1.000 0.707 0.707
## Area1 =~
## X1A 0.200 0.283 0.248
## X1B 0.200 0.283 0.240
## X1C 0.200 0.283 0.296
## X1D 0.200 0.283 0.257
## X1E 0.200 0.283 0.309
## Area2 =~
## X2A 0.125 0.177 0.221
## X2B 0.125 0.177 0.227
## X2C 0.125 0.177 0.227
## X2D 0.125 0.177 0.218
## X2E 0.125 0.177 0.221
## X2F 0.125 0.177 0.219
## X2G 0.125 0.177 0.187
## X2H 0.125 0.177 0.214
## X2_adj 0.205 0.035 5.896 0.000 0.290 0.322
## Area3 =~
## X3A 0.250 0.354 0.375
## X3B 0.250 0.354 0.386
## X3C 0.250 0.354 0.382
## X3D 0.250 0.354 0.363
## Area4 =~
## X4A 0.250 0.354 0.386
## X4B 0.250 0.354 0.472
## X4C 0.250 0.354 0.370
## X4D 0.250 0.354 0.468
## Area5 =~
## X5A 0.250 0.354 0.388
## X5B 0.250 0.354 0.397
## X5C 0.250 0.354 0.480
## X5D 0.250 0.354 0.475
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.000 0.078 0.000 1.000 0.000 0.000
## .X1B -0.009 0.081 -0.111 0.912 -0.009 -0.008
## .X1C -0.001 0.078 -0.009 0.993 -0.001 -0.001
## .X1D 0.000 0.078 0.000 1.000 0.000 0.000
## .X1E 0.001 0.078 0.016 0.987 0.001 0.001
## .X2A -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2C -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2D -0.016 0.082 -0.191 0.849 -0.016 -0.019
## .X2E -0.001 0.078 -0.009 0.993 -0.001 -0.001
## .X2F -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2G -0.003 0.078 -0.044 0.965 -0.003 -0.004
## .X2H -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2_adj 0.000 0.078 0.000 1.000 0.000 0.000
## .X3A -0.001 0.078 -0.016 0.987 -0.001 -0.001
## .X3B 0.000 0.078 0.000 1.000 0.000 0.000
## .X3C 0.000 0.078 0.000 1.000 0.000 0.000
## .X3D 0.000 0.078 0.000 1.000 0.000 0.000
## .X4A 0.000 0.078 0.000 1.000 0.000 0.000
## .X4B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4C 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5A -0.004 0.078 -0.051 0.960 -0.004 -0.004
## .X5B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5C 0.000 0.078 0.000 1.000 0.000 0.000
## .X5D 0.000 0.078 0.000 1.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 1.216 0.117 10.364 0.000 1.216 0.938
## .X1B 1.313 0.137 9.566 0.000 1.313 0.943
## .X1C 0.834 0.100 8.375 0.000 0.834 0.912
## .X1D 1.131 0.110 10.314 0.000 1.131 0.934
## .X1E 0.758 0.096 7.872 0.000 0.758 0.905
## .X2A 0.609 0.065 9.432 0.000 0.609 0.951
## .X2B 0.575 0.067 8.609 0.000 0.575 0.948
## .X2C 0.578 0.067 8.675 0.000 0.578 0.949
## .X2D 0.627 0.058 10.792 0.000 0.627 0.953
## .X2E 0.607 0.071 8.592 0.000 0.607 0.951
## .X2F 0.622 0.070 8.950 0.000 0.622 0.952
## .X2G 0.865 0.108 8.044 0.000 0.865 0.965
## .X2H 0.650 0.065 9.996 0.000 0.650 0.954
## .X2_adj 0.729 0.113 6.463 0.000 0.729 0.896
## .X3A 0.764 0.201 3.799 0.000 0.764 0.859
## .X3B 0.715 0.131 5.465 0.000 0.715 0.851
## .X3C 0.730 0.088 8.281 0.000 0.730 0.854
## .X3D 0.824 0.086 9.581 0.000 0.824 0.868
## .X4A 0.715 0.103 6.936 0.000 0.715 0.851
## .X4B 0.436 0.045 9.585 0.000 0.436 0.777
## .X4C 0.786 0.169 4.645 0.000 0.786 0.863
## .X4D 0.445 0.048 9.277 0.000 0.445 0.781
## .X5A 0.704 0.124 5.693 0.000 0.704 0.849
## .X5B 0.667 0.075 8.840 0.000 0.667 0.842
## .X5C 0.418 0.062 6.720 0.000 0.418 0.770
## .X5D 0.429 0.065 6.611 0.000 0.429 0.774
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.500 0.500
## .Area2 1.000 0.500 0.500
## .Area3 1.000 0.500 0.500
## .Area4 1.000 0.500 0.500
## .Area5 1.000 0.500 0.500# Nested LRT (scaled) — valid because equal is just free + constraints
anova(fit_free, fit_equal)## Warning: lavaan->lavTestLRT():
## some models are based on a different set of observed variables##
## Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
##
## lavaan->lavTestLRT():
## lavaan NOTE: The "Chisq" column contains standard test statistics, not the
## robust test that should be reported per model. A robust difference test is
## a function of two standard (not robust) statistics.
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit_equal 324 11355 11519 2712.4
## fit_free 984 18186 18631 3166.4 444.66 660 1# -------- Plots (pass the FITS, not the strings) ----------
lavaanPlot(model = fit_free, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))lavaanPlot(model = fit_equal, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))Attempt 2
# -------- FREE ----------
mod_free <- '
Economic_Freedom =~ Area1 + Area2 + Area3 + Area4 + Area5
Area2 ~ b_adj*X2_adj
Area1 =~ X1A + X1B + X1C + X1D_i + X1D_ii + X1E
Area2 =~ X2A + X2B + X2C + X2D + X2E + X2F + X2G + X2H
Area3 =~ X3A + X3B + X3C + X3D
Area4 =~ X4A_i + X4A_ii + X4A_iii + X4B_i + X4B_ii + X4C +
X4D_i + X4D_ii + X4D_iii + X4D_iv
Area5 =~ X5A_i + X5A_ii + X5A_iii +
X5B_i + X5B_ii + X5B_iii + X5B_iv + X5B_v + X5B_vi + X5B_vii +
X5C_i + X5C_ii + X5C_iii + X5C_iv +
X5D_i + X5D_ii + X5D_iii
'
# -------- EQUAL ----------
mod_equal <- '
# Top level: each Area contributes equally to Economic_Freedom
Economic_Freedom =~ 1*Area1 + 1*Area2 + 1*Area3 + 1*Area4 + 1*Area5
# Area 2 also has gender adjustment as a predictor
Area2 ~ b_adj*X2_adj
# --- Area 1: 5 indicators → each weight = 1/5 ---
Area1 =~ 0.2*X1A + 0.2*X1B + 0.2*X1C + 0.2*X1D + 0.2*X1E
# --- Area 2: 8 indicators → each weight = 1/8 ---
Area2 =~ 0.125*X2A + 0.125*X2B + 0.125*X2C + 0.125*X2D +
0.125*X2E + 0.125*X2F + 0.125*X2G + 0.125*X2H
# --- Area 3: 4 indicators → each weight = 1/4 ---
Area3 =~ 0.25*X3A + 0.25*X3B + 0.25*X3C + 0.25*X3D
# --- Area 4: 4 indicators → each weight = 1/4 ---
Area4 =~ 0.25*X4A + 0.25*X4B + 0.25*X4C + 0.25*X4D
# --- Area 5: 4 indicators → each weight = 1/4 ---
Area5 =~ 0.25*X5A + 0.25*X5B + 0.25*X5C + 0.25*X5D
'
# -------- Fit & compare ----------
dat <- mydata_std
fit_free <- cfa(mod_free, data = dat, estimator = "MLR", std.lv = TRUE, missing = "fiml")
fit_equal <- cfa(mod_equal, data = dat, estimator = "MLR", std.lv = TRUE, missing = "fiml")
summary(fit_free, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 127 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 141
##
## Number of observations 165
## Number of missing patterns 31
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 3201.363 3041.296
## Degrees of freedom 984 984
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.053
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 6413.070 6013.842
## Degrees of freedom 1035 1035
## P-value 0.000 0.000
## Scaling correction factor 1.066
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.588 0.587
## Tucker-Lewis Index (TLI) 0.566 0.565
##
## Robust Comparative Fit Index (CFI) 0.593
## Robust Tucker-Lewis Index (TLI) 0.572
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -8734.085 -8734.085
## Scaling correction factor 1.198
## for the MLR correction
## Loglikelihood unrestricted model (H1) -7133.404 -7133.404
## Scaling correction factor 1.071
## for the MLR correction
##
## Akaike (AIC) 17750.170 17750.170
## Bayesian (BIC) 18188.108 18188.108
## Sample-size adjusted Bayesian (SABIC) 17741.702 17741.702
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.117 0.113
## 90 Percent confidence interval - lower 0.112 0.108
## 90 Percent confidence interval - upper 0.121 0.117
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.116
## 90 Percent confidence interval - lower 0.111
## 90 Percent confidence interval - upper 0.121
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.126 0.126
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 -1.189 0.342 -3.475 0.001 -0.765 -0.765
## Area2 2.922 0.629 4.647 0.000 0.946 0.946
## Area3 0.450 0.121 3.733 0.000 0.411 0.411
## Area4 2.367 0.493 4.803 0.000 0.921 0.921
## Area5 4.933 2.062 2.392 0.017 0.980 0.980
## Area1 =~
## X1A 0.288 0.062 4.619 0.000 0.447 0.449
## X1B 0.487 0.109 4.475 0.000 0.756 0.761
## X1C -0.296 0.064 -4.651 0.000 -0.459 -0.461
## X1D_i 0.302 0.113 2.667 0.008 0.469 0.470
## X1D_ii 0.274 0.121 2.272 0.023 0.426 0.428
## X1E -0.289 0.066 -4.353 0.000 -0.448 -0.450
## Area2 =~
## X2A 0.292 0.054 5.399 0.000 0.903 0.903
## X2B 0.311 0.055 5.610 0.000 0.959 0.960
## X2C 0.296 0.051 5.822 0.000 0.914 0.914
## X2D 0.249 0.044 5.616 0.000 0.768 0.782
## X2E 0.289 0.051 5.708 0.000 0.891 0.893
## X2F 0.266 0.047 5.630 0.000 0.822 0.823
## X2G 0.099 0.029 3.461 0.001 0.306 0.307
## X2H 0.254 0.044 5.719 0.000 0.783 0.784
## Area3 =~
## X3A 0.647 0.153 4.222 0.000 0.709 0.712
## X3B 0.463 0.110 4.192 0.000 0.508 0.509
## X3C 0.743 0.069 10.735 0.000 0.815 0.818
## X3D 0.109 0.081 1.346 0.178 0.120 0.120
## Area4 =~
## X4A_i 0.220 0.052 4.274 0.000 0.566 0.567
## X4A_ii 0.255 0.056 4.529 0.000 0.656 0.659
## X4A_iii 0.063 0.034 1.830 0.067 0.162 0.163
## X4B_i 0.316 0.061 5.197 0.000 0.813 0.816
## X4B_ii 0.269 0.046 5.890 0.000 0.691 0.693
## X4C 0.143 0.042 3.431 0.001 0.368 0.370
## X4D_i 0.346 0.064 5.419 0.000 0.888 0.891
## X4D_ii 0.197 0.050 3.967 0.000 0.506 0.508
## X4D_iii 0.152 0.038 4.013 0.000 0.391 0.392
## X4D_iv 0.341 0.056 6.114 0.000 0.877 0.879
## Area5 =~
## X5A_i 0.072 0.033 2.189 0.029 0.363 0.364
## X5A_ii 0.039 0.021 1.892 0.059 0.198 0.199
## X5A_iii 0.036 0.023 1.568 0.117 0.183 0.183
## X5B_i 0.080 0.041 1.935 0.053 0.402 0.403
## X5B_ii 0.103 0.050 2.049 0.040 0.520 0.521
## X5B_iii 0.026 0.025 1.029 0.303 0.129 0.129
## X5B_iv 0.019 0.019 0.999 0.318 0.094 0.094
## X5B_v 0.069 0.034 2.043 0.041 0.350 0.350
## X5B_vi 0.026 0.019 1.357 0.175 0.132 0.132
## X5B_vii 0.121 0.053 2.295 0.022 0.611 0.613
## X5C_i 0.079 0.036 2.207 0.027 0.398 0.398
## X5C_ii 0.161 0.063 2.571 0.010 0.810 0.812
## X5C_iii 0.147 0.055 2.697 0.007 0.742 0.747
## X5C_iv 0.137 0.058 2.353 0.019 0.688 0.690
## X5D_i 0.167 0.071 2.359 0.018 0.843 0.846
## X5D_ii 0.041 0.028 1.486 0.137 0.206 0.206
## X5D_iii 0.151 0.067 2.261 0.024 0.758 0.761
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Area2 ~
## X2_adj (b_dj) -0.019 0.134 -0.144 0.885 -0.006 -0.006
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.000 0.078 0.000 1.000 0.000 0.000
## .X1B 0.018 0.078 0.235 0.814 0.018 0.018
## .X1C -0.006 0.079 -0.074 0.941 -0.006 -0.006
## .X1D_i 0.000 0.078 0.000 1.000 0.000 0.000
## .X1D_ii 0.007 0.079 0.086 0.931 0.007 0.007
## .X1E -0.002 0.079 -0.024 0.981 -0.002 -0.002
## .X2A -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2C -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2D -0.063 0.080 -0.783 0.434 -0.063 -0.064
## .X2E -0.003 0.078 -0.036 0.971 -0.003 -0.003
## .X2F -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X2G -0.003 0.078 -0.042 0.966 -0.003 -0.003
## .X2H -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X3A 0.004 0.077 0.048 0.961 0.004 0.004
## .X3B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X3C -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X3D -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4A_i -0.024 0.081 -0.298 0.765 -0.024 -0.024
## .X4A_ii -0.001 0.078 -0.018 0.986 -0.001 -0.001
## .X4A_iii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4B_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4B_ii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4C -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4D_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4D_ii 0.003 0.078 0.039 0.969 0.003 0.003
## .X4D_iii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X4D_iv -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5A_i -0.029 0.084 -0.342 0.733 -0.029 -0.029
## .X5A_ii -0.005 0.079 -0.061 0.952 -0.005 -0.005
## .X5A_iii -0.004 0.079 -0.047 0.962 -0.004 -0.004
## .X5B_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5B_ii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5B_iii -0.017 0.080 -0.212 0.832 -0.017 -0.017
## .X5B_iv -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5B_v -0.007 0.079 -0.095 0.925 -0.007 -0.007
## .X5B_vi -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5B_vii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5C_i -0.019 0.079 -0.246 0.806 -0.019 -0.019
## .X5C_ii -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5C_iii 0.003 0.078 0.037 0.971 0.003 0.003
## .X5C_iv -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5D_i -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5D_ii -0.008 0.079 -0.107 0.915 -0.008 -0.008
## .X5D_iii -0.000 0.078 -0.000 1.000 -0.000 -0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.794 0.126 6.301 0.000 0.794 0.799
## .X1B 0.415 0.099 4.187 0.000 0.415 0.420
## .X1C 0.781 0.089 8.811 0.000 0.781 0.787
## .X1D_i 0.774 0.101 7.650 0.000 0.774 0.779
## .X1D_ii 0.812 0.097 8.402 0.000 0.812 0.817
## .X1E 0.791 0.094 8.374 0.000 0.791 0.797
## .X2A 0.184 0.025 7.317 0.000 0.184 0.184
## .X2B 0.078 0.017 4.548 0.000 0.078 0.079
## .X2C 0.164 0.033 4.995 0.000 0.164 0.164
## .X2D 0.374 0.044 8.575 0.000 0.374 0.388
## .X2E 0.201 0.027 7.520 0.000 0.201 0.202
## .X2F 0.321 0.042 7.671 0.000 0.321 0.322
## .X2G 0.901 0.107 8.410 0.000 0.901 0.906
## .X2H 0.384 0.053 7.231 0.000 0.384 0.385
## .X3A 0.488 0.138 3.529 0.000 0.488 0.492
## .X3B 0.736 0.149 4.936 0.000 0.736 0.740
## .X3C 0.329 0.096 3.416 0.001 0.329 0.331
## .X3D 0.980 0.091 10.736 0.000 0.980 0.985
## .X4A_i 0.675 0.121 5.578 0.000 0.675 0.678
## .X4A_ii 0.561 0.117 4.810 0.000 0.561 0.566
## .X4A_iii 0.968 0.154 6.300 0.000 0.968 0.974
## .X4B_i 0.333 0.058 5.775 0.000 0.333 0.335
## .X4B_ii 0.516 0.068 7.549 0.000 0.516 0.519
## .X4C 0.858 0.186 4.611 0.000 0.858 0.863
## .X4D_i 0.206 0.041 5.073 0.000 0.206 0.207
## .X4D_ii 0.738 0.103 7.180 0.000 0.738 0.742
## .X4D_iii 0.841 0.084 9.982 0.000 0.841 0.846
## .X4D_iv 0.226 0.040 5.576 0.000 0.226 0.227
## .X5A_i 0.865 0.093 9.319 0.000 0.865 0.868
## .X5A_ii 0.955 0.138 6.923 0.000 0.955 0.961
## .X5A_iii 0.961 0.152 6.316 0.000 0.961 0.966
## .X5B_i 0.832 0.086 9.661 0.000 0.832 0.838
## .X5B_ii 0.724 0.082 8.861 0.000 0.724 0.728
## .X5B_iii 0.977 0.178 5.480 0.000 0.977 0.983
## .X5B_iv 0.985 0.106 9.277 0.000 0.985 0.991
## .X5B_v 0.876 0.089 9.802 0.000 0.876 0.877
## .X5B_vi 0.977 0.056 17.379 0.000 0.977 0.982
## .X5B_vii 0.621 0.086 7.242 0.000 0.621 0.624
## .X5C_i 0.840 0.093 9.036 0.000 0.840 0.841
## .X5C_ii 0.338 0.039 8.563 0.000 0.338 0.340
## .X5C_iii 0.436 0.074 5.905 0.000 0.436 0.442
## .X5C_iv 0.520 0.068 7.643 0.000 0.520 0.524
## .X5D_i 0.283 0.042 6.687 0.000 0.283 0.285
## .X5D_ii 0.953 0.127 7.495 0.000 0.953 0.957
## .X5D_iii 0.419 0.071 5.927 0.000 0.419 0.422
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.414 0.414
## .Area2 1.000 0.105 0.105
## .Area3 1.000 0.831 0.831
## .Area4 1.000 0.151 0.151
## .Area5 1.000 0.039 0.039summary(fit_equal, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 20 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 51
##
## Number of observations 165
## Number of missing patterns 13
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2608.338 2362.185
## Degrees of freedom 324 324
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.104
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 3473.006 3108.849
## Degrees of freedom 325 325
## P-value 0.000 0.000
## Scaling correction factor 1.117
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.274 0.268
## Tucker-Lewis Index (TLI) 0.272 0.266
##
## Robust Comparative Fit Index (CFI) 0.278
## Robust Tucker-Lewis Index (TLI) 0.275
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5338.721 -5338.721
## Scaling correction factor 1.258
## for the MLR correction
## Loglikelihood unrestricted model (H1) -4034.552 -4034.552
## Scaling correction factor 1.125
## for the MLR correction
##
## Akaike (AIC) 10779.441 10779.441
## Bayesian (BIC) 10937.844 10937.844
## Sample-size adjusted Bayesian (SABIC) 10776.378 10776.378
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.207 0.195
## 90 Percent confidence interval - lower 0.199 0.188
## 90 Percent confidence interval - upper 0.214 0.202
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.206
## 90 Percent confidence interval - lower 0.198
## 90 Percent confidence interval - upper 0.214
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.352 0.352
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Economic_Freedom =~
## Area1 1.000 0.707 0.707
## Area2 1.000 0.362 0.362
## Area3 1.000 0.707 0.707
## Area4 1.000 0.707 0.707
## Area5 1.000 0.707 0.707
## Area1 =~
## X1A 0.200 0.283 0.249
## X1B 0.200 0.283 0.242
## X1C 0.200 0.283 0.295
## X1D 0.200 0.283 0.259
## X1E 0.200 0.283 0.305
## Area2 =~
## X2A 0.125 0.346 0.430
## X2B 0.125 0.346 0.446
## X2C 0.125 0.346 0.443
## X2D 0.125 0.346 0.431
## X2E 0.125 0.346 0.435
## X2F 0.125 0.346 0.421
## X2G 0.125 0.346 0.344
## X2H 0.125 0.346 0.401
## Area3 =~
## X3A 0.250 0.354 0.380
## X3B 0.250 0.354 0.382
## X3C 0.250 0.354 0.387
## X3D 0.250 0.354 0.362
## Area4 =~
## X4A 0.250 0.354 0.382
## X4B 0.250 0.354 0.464
## X4C 0.250 0.354 0.371
## X4D 0.250 0.354 0.463
## Area5 =~
## X5A 0.250 0.354 0.384
## X5B 0.250 0.354 0.396
## X5C 0.250 0.354 0.477
## X5D 0.250 0.354 0.471
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Area2 ~
## X2_adj (b_dj) 2.384 0.425 5.609 0.000 0.862 0.859
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 0.000 0.078 0.000 1.000 0.000 0.000
## .X1B -0.008 0.081 -0.094 0.925 -0.008 -0.007
## .X1C 0.000 0.078 0.002 0.999 0.000 0.000
## .X1D 0.000 0.078 0.000 1.000 0.000 0.000
## .X1E 0.001 0.078 0.018 0.986 0.001 0.002
## .X2A 0.000 0.072 0.000 1.000 0.000 0.000
## .X2B 0.000 0.072 0.000 1.000 0.000 0.000
## .X2C 0.000 0.072 0.000 1.000 0.000 0.000
## .X2D -0.019 0.075 -0.259 0.796 -0.019 -0.024
## .X2E -0.001 0.071 -0.010 0.992 -0.001 -0.001
## .X2F 0.000 0.073 0.000 1.000 0.000 0.000
## .X2G -0.004 0.077 -0.052 0.958 -0.004 -0.004
## .X2H -0.000 0.075 -0.000 1.000 -0.000 -0.000
## .X3A -0.001 0.078 -0.007 0.994 -0.001 -0.001
## .X3B 0.000 0.078 0.000 1.000 0.000 0.000
## .X3C 0.000 0.078 0.000 1.000 0.000 0.000
## .X3D 0.000 0.078 0.000 1.000 0.000 0.000
## .X4A 0.000 0.078 0.000 1.000 0.000 0.000
## .X4B 0.000 0.078 0.000 1.000 0.000 0.000
## .X4C 0.000 0.078 0.000 1.000 0.000 0.000
## .X4D 0.000 0.078 0.000 1.000 0.000 0.000
## .X5A -0.004 0.078 -0.047 0.963 -0.004 -0.004
## .X5B -0.000 0.078 -0.000 1.000 -0.000 -0.000
## .X5C 0.000 0.078 0.000 1.000 0.000 0.000
## .X5D -0.000 0.078 -0.000 1.000 -0.000 -0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1A 1.214 0.118 10.297 0.000 1.214 0.938
## .X1B 1.285 0.135 9.535 0.000 1.285 0.941
## .X1C 0.841 0.100 8.391 0.000 0.841 0.913
## .X1D 1.111 0.108 10.336 0.000 1.111 0.933
## .X1E 0.782 0.097 8.059 0.000 0.782 0.907
## .X2A 0.527 0.056 9.490 0.000 0.527 0.815
## .X2B 0.480 0.054 8.949 0.000 0.480 0.801
## .X2C 0.489 0.057 8.652 0.000 0.489 0.804
## .X2D 0.525 0.051 10.396 0.000 0.525 0.815
## .X2E 0.511 0.063 8.086 0.000 0.511 0.810
## .X2F 0.555 0.059 9.345 0.000 0.555 0.823
## .X2G 0.888 0.110 8.054 0.000 0.888 0.881
## .X2H 0.623 0.065 9.606 0.000 0.623 0.839
## .X3A 0.741 0.204 3.639 0.000 0.741 0.856
## .X3B 0.731 0.133 5.494 0.000 0.731 0.854
## .X3C 0.710 0.088 8.074 0.000 0.710 0.850
## .X3D 0.829 0.087 9.517 0.000 0.829 0.869
## .X4A 0.731 0.104 7.059 0.000 0.731 0.854
## .X4B 0.455 0.047 9.611 0.000 0.455 0.785
## .X4C 0.783 0.173 4.536 0.000 0.783 0.862
## .X4D 0.459 0.049 9.383 0.000 0.459 0.786
## .X5A 0.724 0.130 5.558 0.000 0.724 0.853
## .X5B 0.672 0.076 8.886 0.000 0.672 0.843
## .X5C 0.424 0.064 6.642 0.000 0.424 0.772
## .X5D 0.438 0.065 6.698 0.000 0.438 0.778
## Economic_Fredm 1.000 1.000 1.000
## .Area1 1.000 0.500 0.500
## .Area2 1.000 0.131 0.131
## .Area3 1.000 0.500 0.500
## .Area4 1.000 0.500 0.500
## .Area5 1.000 0.500 0.500# Nested LRT (scaled) — valid because equal is just free + constraints
anova(fit_free, fit_equal)## Warning: lavaan->lavTestLRT():
## some models are based on a different set of observed variables##
## Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
##
## lavaan->lavTestLRT():
## lavaan NOTE: The "Chisq" column contains standard test statistics, not the
## robust test that should be reported per model. A robust difference test is
## a function of two standard (not robust) statistics.
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit_equal 324 10779 10938 2608.3
## fit_free 984 17750 18188 3201.4 577.26 660 0.9909# -------- Plots (pass the FITS, not the strings) ----------
lavaanPlot(model = fit_free, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))lavaanPlot(model = fit_equal, coefs = TRUE, covs = FALSE,
graph_options = list(rankdir = "LR"),
edge_options = list(color = "grey40"))EFA
set.seed(42) # for stable parallel analysis plots
# --- 0) Helper: pretty names from stored labels ------------------------------
pretty_names <- function(df, vars) {
sapply(vars, function(v) {
lb <- attr(df[[v]], "label", exact = TRUE)
if (is.null(lb) || is.na(lb) || lb == "") v else lb
}, USE.NAMES = FALSE)
}
# --- 1) Build the EFA data matrix (numeric only, selected columns) -----------
efa_cols <- c(
"X1A","X1B","X1C","X1D_i","X1D_ii","X1E",
"X2A","X2B","X2C","X2D","X2E","X2F","X2G","X2H",
"X3A","X3B","X3C","X3D",
"X4A_i","X4A_ii","X4A_iii","X4B_i","X4B_ii","X4C",
"X4D_i","X4D_ii","X4D_iii","X4D_iv",
"X5A_i","X5A_ii","X5A_iii",
"X5B_i","X5B_ii","X5B_iii","X5B_iv","X5B_v","X5B_vi","X5B_vii",
"X5C_i","X5C_ii","X5C_iii","X5C_iv",
"X5D_i","X5D_ii","X5D_iii"
)
# keep only columns that exist and are numeric
efa_vars <- mydata_std %>%
dplyr::select(dplyr::any_of(efa_cols)) %>%
dplyr::select(where(is.numeric))
# labels aligned to the actual efa_vars we ended up with
var_labels <- pretty_names(mydata_std, colnames(efa_vars))
# --- 2) Factorability checks -------------------------------------------------
# KMO (≥ .60 is ok)
print(psych::KMO(efa_vars))## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = efa_vars)
## Overall MSA = 0.86
## MSA for each item =
## X1A X1B X1C X1D_i X1D_ii X1E X2A X2B X2C X2D
## 0.84 0.91 0.81 0.70 0.59 0.83 0.91 0.91 0.93 0.89
## X2E X2F X2G X2H X3A X3B X3C X3D X4A_i X4A_ii
## 0.89 0.95 0.71 0.92 0.70 0.68 0.69 0.88 0.85 0.86
## X4A_iii X4B_i X4B_ii X4C X4D_i X4D_ii X4D_iii X4D_iv X5A_i X5A_ii
## 0.46 0.93 0.92 0.80 0.92 0.78 0.85 0.94 0.70 0.78
## X5A_iii X5B_i X5B_ii X5B_iii X5B_iv X5B_v X5B_vi X5B_vii X5C_i X5C_ii
## 0.57 0.71 0.74 0.74 0.39 0.59 0.41 0.89 0.81 0.93
## X5C_iii X5C_iv X5D_i X5D_ii X5D_iii
## 0.93 0.93 0.92 0.63 0.90# Bartlett (p < .05 => correlations not identity)
R <- cor(efa_vars, use = "pairwise.complete.obs")
print(psych::cortest.bartlett(R, n = nrow(efa_vars)))## $chisq
## [1] 5867.33
##
## $p.value
## [1] 0
##
## $df
## [1] 990# --- 3) Parallel analysis (suggested # of factors) ---------------------------
fa.parallel(efa_vars, fm = "ml", fa = "fa", use = "pairwise")
## Parallel analysis suggests that the number of factors = 7 and the number of components = NA# --- 4) EFA: 5 factors, ML, oblimin -----------------------------------------
efa_result <- fa(efa_vars, nfactors = 5, fm = "ml", rotate = "oblimin",
use = "pairwise")
# quick console print (cutoff .30)
print(efa_result$loadings, cutoff = 0.30)##
## Loadings:
## ML2 ML3 ML5 ML1 ML4
## X1A -0.600 0.545
## X1B -0.402 -0.393 0.438
## X1C
## X1D_i -0.563 0.301
## X1D_ii 0.531
## X1E 0.635
## X2A 0.870
## X2B 0.847
## X2C 0.608 0.328
## X2D 0.717
## X2E 0.893
## X2F 0.538
## X2G 0.483
## X2H 0.601
## X3A 0.331
## X3B 0.470
## X3C 0.320
## X3D 0.670
## X4A_i 0.744
## X4A_ii 0.644
## X4A_iii 0.357
## X4B_i 0.665
## X4B_ii 0.372 0.385
## X4C
## X4D_i 0.630
## X4D_ii 0.471
## X4D_iii
## X4D_iv 0.313 0.446 0.417
## X5A_i 0.328
## X5A_ii
## X5A_iii
## X5B_i 0.468
## X5B_ii 0.980
## X5B_iii -0.344 0.391
## X5B_iv
## X5B_v 0.914
## X5B_vi
## X5B_vii 0.373 0.400
## X5C_i 0.729
## X5C_ii 0.555
## X5C_iii 0.810
## X5C_iv 0.349 0.378
## X5D_i 0.391 0.376
## X5D_ii 0.479
## X5D_iii 0.313 0.597
##
## ML2 ML3 ML5 ML1 ML4
## SS loadings 7.016 4.530 2.564 2.283 2.527
## Proportion Var 0.156 0.101 0.057 0.051 0.056
## Cumulative Var 0.156 0.257 0.314 0.364 0.420# --- 5) Clean loading table with labels --------------------------------------
load_tbl <- as.data.frame(unclass(efa_result$loadings))
load_tbl <- tibble::rownames_to_column(load_tbl, "var")
colnames(load_tbl) <- c("var", paste0("F", seq_len(ncol(load_tbl) - 1L)))
# attach labels
load_tbl$label <- var_labels[match(load_tbl$var, colnames(efa_vars))]
# primary factor (largest |loading|)
Fcols <- paste0("F", 1:5)
load_tbl$primary <- apply(load_tbl[, Fcols, drop = FALSE], 1, function(x) {
paste0("F", which.max(abs(x)))
})
load_tbl$loading <- apply(load_tbl[, Fcols, drop = FALSE], 1, function(x) {
x[which.max(abs(x))]
})
cat("\n\n=== Primary loadings (|λ| ≥ .30) ===\n")##
##
## === Primary loadings (|λ| ≥ .30) ===print(
load_tbl %>%
mutate(keep = abs(loading) >= .30) %>%
filter(keep) %>%
arrange(primary, desc(abs(loading))) %>%
select(var, label, primary, loading),
digits = 3, row.names = FALSE
)## var label primary loading
## X2E X2E F1 0.893
## X2A X2A F1 0.870
## X2B X2B F1 0.847
## X5C_iii X5C_iii F1 0.810
## X2D X2D F1 0.717
## X2C X2C F1 0.608
## X2H X2H F1 0.601
## X1A X1A F1 -0.600
## X1D_i X1D_i F1 -0.563
## X5C_ii X5C_ii F1 0.555
## X2F X2F F1 0.538
## X5D_i X5D_i F1 0.391
## X3C X3C F1 0.320
## X4A_i X4A_i F2 0.744
## X3D X3D F2 0.670
## X4B_i X4B_i F2 0.665
## X4A_ii X4A_ii F2 0.644
## X4D_i X4D_i F2 0.630
## X2G X2G F2 0.483
## X4D_ii X4D_ii F2 0.471
## X4D_iv X4D_iv F2 0.446
## X4B_ii X4B_ii F2 0.385
## X4A_iii X4A_iii F2 0.357
## X1E X1E F3 0.635
## X5D_iii X5D_iii F3 0.597
## X3B X3B F3 0.470
## X5B_vii X5B_vii F3 0.400
## X5A_i X5A_i F3 0.328
## X5B_ii X5B_ii F4 0.980
## X5B_v X5B_v F4 0.914
## X5C_i X5C_i F5 0.729
## X1D_ii X1D_ii F5 0.531
## X5D_ii X5D_ii F5 0.479
## X5B_i X5B_i F5 0.468
## X1B X1B F5 0.438
## X5B_iii X5B_iii F5 0.391
## X5C_iv X5C_iv F5 0.378
## X3A X3A F5 0.331# --- 6) Diagram (psych) ------------------------------------------------------
#fa.diagram(efa_result, main = "Exploratory Factor Analysis (5 Factors)")
# --- 7) Heatmap of loadings (clean, grouped) ---------------------------------
# tidy long
long_load <- load_tbl %>%
select(var, label, all_of(Fcols)) %>%
tidyr::pivot_longer(all_of(Fcols), names_to = "Factor", values_to = "Loading")
# create a stable item order: by primary factor, then by |loading|
item_order <- load_tbl %>%
arrange(primary, desc(abs(loading))) %>%
pull(label)
long_load$label <- factor(long_load$label, levels = rev(unique(item_order))) # top-to-bottom
ggplot(long_load, aes(x = Factor, y = label, fill = Loading)) +
geom_tile() +
scale_fill_gradient2(limits = c(-1, 1)) +
labs(x = NULL, y = NULL, title = "Factor Loadings Heatmap (ML, oblimin)") +
theme_minimal(base_size = 11) +
theme(
panel.grid = element_blank(),
axis.text.y = element_text(size = 8)
)
# --- 8) (Optional) Factor scores for later use --------------------------------
scores <- factor.scores(efa_vars, efa_result, method = "tenBerge")$scores
head(scores)## ML2 ML3 ML5 ML1 ML4
## [1,] -0.2334360 0.5775833 0.83542056 -0.08241388 -0.1291864
## [2,] -0.6439542 -0.9792912 -1.89988411 0.22915325 0.1604326
## [3,] -1.0835857 -0.8704263 -1.28908081 -0.46166276 -0.2248049
## [4,] 0.2475448 -0.8425737 -0.55756147 -1.74652397 -2.2407401
## [5,] -0.2699180 0.3825257 0.07215822 0.68381881 0.5855438
## [6,] 1.9081423 0.7283472 0.45109573 0.84061607 0.1717790
# Did both models converge?
lavaan::inspect(fit_free, "converged")## [1] TRUElavaan::inspect(fit_equal, "converged")## [1] TRUE# Degrees of freedom and robust scaling factors available?
fitMeasures(fit_free, c("df","chisq","pvalue","scaling.factor.h0"))## df chisq pvalue scaling.factor.h0
## 984.000 3201.363 0.000 1.198fitMeasures(fit_equal, c("df","chisq","pvalue","scaling.factor.h0"))## df chisq pvalue scaling.factor.h0
## 324.000 2608.338 0.000 1.258anova(fit_free, fit_equal)## Warning: lavaan->lavTestLRT():
## some models are based on a different set of observed variables##
## Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
##
## lavaan->lavTestLRT():
## lavaan NOTE: The "Chisq" column contains standard test statistics, not the
## robust test that should be reported per model. A robust difference test is
## a function of two standard (not robust) statistics.
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit_equal 324 10779 10938 2608.3
## fit_free 984 17750 18188 3201.4 577.26 660 0.9909lavaan::lavTestLRT(fit_free, fit_equal, method = "satorra.bentler.2010")## Warning: lavaan->lavTestLRT():
## some models are based on a different set of observed variables##
## Scaled Chi-Squared Difference Test (method = "satorra.bentler.2010")
##
## lavaan->lavTestLRT():
## lavaan NOTE: The "Chisq" column contains standard test statistics, not the
## robust test that should be reported per model. A robust difference test is
## a function of two standard (not robust) statistics.
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit_equal 324 10779 10938 2608.3
## fit_free 984 17750 18188 3201.4 550.93 660 0.9992# alternatives you can try if 2010 still fails:
# lavaan::lavTestLRT(fit_free, fit_equal, method = "satorra.bentler")
# lavaan::lavTestLRT(fit_free, fit_equal, method = "yuan.bentler")
# lavaan::lavTestLRT(fit_free, fit_equal, method = "scaled.shifted")fit_free_ML <- cfa(mod_free, data = mydata_std, estimator = "ML", std.lv = TRUE, missing = "fiml")
fit_equal_ML <- cfa(mod_equal, data = mydata_std, estimator = "ML", std.lv = TRUE, missing = "fiml")
anova(fit_free_ML, fit_equal_ML)## Warning: lavaan->lavTestLRT():
## some models are based on a different set of observed variables##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_equal_ML 324 10779 10938 2608.3
## fit_free_ML 984 17750 18188 3201.4 593.03 0 660 0.9707Citations
citation("lavaan")## To cite lavaan in publications use:
##
## Yves Rosseel (2012). lavaan: An R Package for Structural Equation
## Modeling. Journal of Statistical Software, 48(2), 1-36.
## https://doi.org/10.18637/jss.v048.i02
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {{lavaan}: An {R} Package for Structural Equation Modeling},
## author = {Yves Rosseel},
## journal = {Journal of Statistical Software},
## year = {2012},
## volume = {48},
## number = {2},
## pages = {1--36},
## doi = {10.18637/jss.v048.i02},
## }