CIFA Urban Ipixuna: model fitting
Perform confirmatory item factor analysis (CIFA) using the results obtained from the exploratory item factor analysis (EIFA) with 3 latent dimensions.
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, unionlibrary(mirt)#> Loading required package: stats4#> Loading required package: latticelibrary(spifa)
library(sf)#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUElibrary(purrr)
library(ggplot2)
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")
fidata <- file.path(path_processed, "fi-items-ipixuna-urban.gpkg") |>
st_read(as_tibble = TRUE)#> Reading layer `fi-items-ipixuna-urban' from data source
#> `/home/rstudio/documents/projects/food-insecurity-mapping/data/processed/fi-items-ipixuna-urban.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 200 features and 36 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -71.70038 ymin: -7.06058 xmax: -71.68109 ymax: -7.03724
#> Geodetic CRS: WGS 84eifa_data <- readRDS(file.path(path_modelled, "eifa-ipixuna-urban.rds"))Get items data and eifa model
# items
items_data <- fidata |>
st_set_geometry(NULL) |>
dplyr::select(matches("^item_[0-9]+_[A-D]"))
# discrimination parameters with rotation
eifa_discrimiation <- summary(eifa_data$model[[3]], rotate = "varimax", suppres = 0.5)#>
#> Rotation: varimax
#>
#> Rotated factor loadings:
#>
#> F1 F2 F3 h2
#> item_1_A_worried_that_food_ends NA -0.639 NA 0.742
#> item_2_A_run_out_of_food NA NA -0.529 0.689
#> item_3_A_ate_few_food_types -0.761 NA NA 0.831
#> item_4_B_skipped_a_meal -0.764 NA -0.506 0.906
#> item_5_B_ate_less_than_required -0.630 NA -0.653 0.910
#> item_6_B_hungry_but_did_not_eat -0.683 NA -0.634 0.993
#> item_7_B_at_most_one_meal_per_day -0.736 NA NA 0.817
#> item_8_C_ate_few_food_types -0.836 NA NA 0.825
#> item_9_C_ate_less_than_required -0.858 NA NA 0.880
#> item_10_C_decreased_food_quantity -0.875 NA NA 0.926
#> item_11_C_skipped_a_meal -0.935 NA NA 0.985
#> item_12_C_hungry_but_did_not_eat -0.881 NA NA 0.949
#> item_13_C_at_most_one_meal_per_day -0.946 NA NA 0.994
#> item_14_D_food_just_with_farinha -0.568 NA -0.640 0.797
#> item_15_D_credit_for_eating NA -0.753 NA 0.668
#> item_16_D_borrowed_food NA -0.892 NA 0.990
#> item_17_D_had_meals_at_neighbors NA -0.552 NA 0.581
#> item_18_D_reduced_meat_or_fish -0.541 NA NA 0.634
#>
#> Rotated SS loadings: 8.571 3.31 3.235
#>
#> Factor correlations:
#>
#> F1 F2 F3
#> F1 1 0 0
#> F2 0 1 0
#> F3 0 0 1Execute confirmatory item factor analysis
# general dimensions
q <- ncol(items_data) # number of items
m <- 3 # number of latent factors
# restrictions on the discrimination parameters A for confirmatory analysis
L_a <- (abs(eifa_discrimiation$rotF) > 0.5) * 1
# prior for the discrimination parameters A
A_mean <- matrix(0, q, m)
A_mean[c(11, 13), 1] <- 1
A_mean[c(16), 2] <- 1
A_mean[c(14), 3] <- 1
A_sd <- matrix(1, q, m)
A_sd[A_mean == 1] <- 0.45
# execute model
iter <- 1.5 * 10 ^ 6
thin <- 200
samples <- spifa(
response = item_1_A_worried_that_food_ends:item_18_D_reduced_meat_or_fish,
data = items_data, nfactors = 3, niter = iter, thin = thin,
constrains = list(A = L_a),
A_opt = list(initial = A_mean, prior_mean = A_mean, prior_sd = A_sd)
)Print the settings defined for the model.
attr(samples, "model_info")[-c(1:4)]#> $nobs
#> [1] 200
#>
#> $nitems
#> [1] 18
#>
#> $nfactors
#> [1] 3
#>
#> $ngp
#> [1] 3
#>
#> $niter
#> [1] 1500000
#>
#> $thin
#> [1] 200
#>
#> $standardize
#> [1] TRUE
#>
#> $constrain_L
#> F1 F2 F3
#> item_1_A_worried_that_food_ends 0 1 0
#> item_2_A_run_out_of_food 0 0 1
#> item_3_A_ate_few_food_types 1 0 0
#> item_4_B_skipped_a_meal 1 0 1
#> item_5_B_ate_less_than_required 1 0 1
#> item_6_B_hungry_but_did_not_eat 1 0 1
#> item_7_B_at_most_one_meal_per_day 1 0 0
#> item_8_C_ate_few_food_types 1 0 0
#> item_9_C_ate_less_than_required 1 0 0
#> item_10_C_decreased_food_quantity 1 0 0
#> item_11_C_skipped_a_meal 1 0 0
#> item_12_C_hungry_but_did_not_eat 1 0 0
#> item_13_C_at_most_one_meal_per_day 1 0 0
#> item_14_D_food_just_with_farinha 1 0 1
#> item_15_D_credit_for_eating 0 1 0
#> item_16_D_borrowed_food 0 1 0
#> item_17_D_had_meals_at_neighbors 0 1 0
#> item_18_D_reduced_meat_or_fish 1 0 0
#>
#> $constrain_T
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
#>
#> $constrain_V_sd
#> [,1]
#> [1,] 1
#> [2,] 1
#> [3,] 1
#>
#> $adap_Sigma
#> [,1] [,2] [,3]
#> [1,] 0.001 0.000 0.000
#> [2,] 0.000 0.001 0.000
#> [3,] 0.000 0.000 0.001
#>
#> $adap_scale
#> [1] 1
#>
#> $adap_C
#> [1] 0.7
#>
#> $adap_alpha
#> [1] 0.8
#>
#> $adap_accep_prob
#> [1] 0.234
#>
#> $c_initial
#> [1] -0.458832737 -1.598318595 0.002298863 -0.114153014 1.257591637 -0.952911389 1.137376323
#> [8] -1.260219269 0.780287540 -0.247557977 -0.844603367 -0.699546482 0.146197106 -0.816330230
#> [15] 0.673824769 -0.290361446 0.476389343 -0.157801694
#>
#> $c_prior_mean
#> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>
#> $c_prior_sd
#> [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#>
#> $A_initial
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 0 0 0
#> [3,] 0 0 0
#> [4,] 0 0 0
#> [5,] 0 0 0
#> [6,] 0 0 0
#> [7,] 0 0 0
#> [8,] 0 0 0
#> [9,] 0 0 0
#> [10,] 0 0 0
#> [11,] 1 0 0
#> [12,] 0 0 0
#> [13,] 1 0 0
#> [14,] 0 0 1
#> [15,] 0 0 0
#> [16,] 0 1 0
#> [17,] 0 0 0
#> [18,] 0 0 0
#>
#> $A_prior_mean
#> [,1] [,2] [,3]
#> [1,] 0 0 0
#> [2,] 0 0 0
#> [3,] 0 0 0
#> [4,] 0 0 0
#> [5,] 0 0 0
#> [6,] 0 0 0
#> [7,] 0 0 0
#> [8,] 0 0 0
#> [9,] 0 0 0
#> [10,] 0 0 0
#> [11,] 1 0 0
#> [12,] 0 0 0
#> [13,] 1 0 0
#> [14,] 0 0 1
#> [15,] 0 0 0
#> [16,] 0 1 0
#> [17,] 0 0 0
#> [18,] 0 0 0
#>
#> $A_prior_sd
#> [,1] [,2] [,3]
#> [1,] 1.00 1.00 1.00
#> [2,] 1.00 1.00 1.00
#> [3,] 1.00 1.00 1.00
#> [4,] 1.00 1.00 1.00
#> [5,] 1.00 1.00 1.00
#> [6,] 1.00 1.00 1.00
#> [7,] 1.00 1.00 1.00
#> [8,] 1.00 1.00 1.00
#> [9,] 1.00 1.00 1.00
#> [10,] 1.00 1.00 1.00
#> [11,] 0.45 1.00 1.00
#> [12,] 1.00 1.00 1.00
#> [13,] 0.45 1.00 1.00
#> [14,] 1.00 1.00 0.45
#> [15,] 1.00 1.00 1.00
#> [16,] 1.00 0.45 1.00
#> [17,] 1.00 1.00 1.00
#> [18,] 1.00 1.00 1.00
#>
#> $R_initial
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
#>
#> $R_prior_eta
#> [1] 1.5
#>
#> $B_initial
#> [,1] [,2] [,3]
#> [1,] NA NA NA
#>
#> $B_prior_mean
#> [,1] [,2] [,3]
#> [1,] NA NA NA
#>
#> $B_prior_sd
#> [,1] [,2] [,3]
#> [1,] NA NA NA
#>
#> $sigmas_gp_initial
#> [1] NA NA NA
#>
#> $sigmas_gp_mean
#> [1] NA NA NA
#>
#> $sigmas_gp_sd
#> [1] NA NA NA
#>
#> $phi_gp_initial
#> [1] NA NA NA
#>
#> $phi_gp_mean
#> [1] NA NA NA
#>
#> $phi_gp_sd
#> [1] NA NA NA
#>
#> $model_type
#> [1] "cifa"Save model output
saveRDS(samples, file = file.path(path_modelled, "cifa-ipixuna-urban.rds"))Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 5685.509 9907.787 1975.831