EIFA Urban Ipixuna: compare models

Compare EIFA models created with different numbers of dimension, so we can select the number of dimensions for the confirmatory factor analysis and the spatial factor analysis.

Load required libraries and data

rm(list = ls())
library(day2day)
library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(mirt)

#> Loading required package: stats4

#> Loading required package: lattice

library(sf)

#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(purrr)

path_main <- git_path()
path_data <- file.path(path_main, "data")
path_raw <- file.path(path_data, "raw")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")

eifa_data <- readRDS(file.path(path_modelled, "eifa-ipixuna-urban.rds"))

Compare models with anova

Obtain table comparing all models.

ndims <- 2:nrow(eifa_data)
test_table <- map(ndims, ~anova(eifa_data$model[[.-1]], eifa_data$model[[.]])) |>
    reduce(rbind) |>
    slice(c(1, 1:length(ndims) * 2)) |>
    mutate(across(AIC:X2, ~ round(.x, 2))) |>
    mutate(p = round(p, 3)) |>
    mutate(test = c(NaN, paste0(ndims - 1, " vs ", ndims))) |>
    mutate(model = 1:nrow(eifa_data)) |>
    relocate(model)

#> 
#> Model 1: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> Model 2: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> 
#> 
#> Model 1: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> Model 2: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> 
#> 
#> Model 1: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> Model 2: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> 
#> 
#> Model 1: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> Model 2: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> 
#> 
#> Model 1: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))
#> Model 2: mirt(data = items_data, model = x, method = "MHRM", verbose = FALSE, 
#>     technical = list(NCYCLES = 20000))

Visualize table of model comparison

knitr::kable(test_table)

model	AIC	AICc	SABIC	HQ	BIC	logLik	X2	df	p	test
1	2650.77	2667.11	2655.45	2698.82	2769.51	-1289.38	NaN	NaN	NaN	NaN
2	2614.86	2654.07	2621.76	2685.60	2789.67	-1254.43	69.91	17	0.000	1 vs 2
3	2621.42	2695.73	2630.41	2713.52	2849.01	-1241.71	25.44	16	0.062	2 vs 3
4	2638.32	2762.49	2649.26	2750.44	2915.38	-1235.16	13.10	15	0.594	3 vs 4
5	2645.92	2838.03	2658.68	2776.72	2969.15	-1224.96	20.41	14	0.118	4 vs 5
6	2694.61	2977.15	2709.06	2842.77	3060.72	-1236.30	-22.69	13	1.000	5 vs 6

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   1.671   0.348   1.700