Integrated workflow
vignette("spifa")
fits the full spatial model directly, using a discrimination restriction
and a spatial range prior already decided. In practice both could be
informed by previous knowledge or a preliminary analysis on the data
itself. This vignette walks through that preliminary analysis, then fits
the full model:
- Perform non-spatial exploratory item factor analysis to suggest a discrimination structure using the mirt package.
- Analyse the empirical variogram of that fit’s factor scores to suggest a spatial range using the gstat package.
- Fit the full spatial model with
spifa(), using the restriction and range prior from the previous two steps.
#> Loading required package: stats4#> Loading required package: latticelibrary(gstat)
library(sf)
library(spifa)
library(ggplot2)
data(ipixuna)
nitems <- ncol(ipixuna$items)
nfactors <- 3Step 1: Discrimination structure with mirt
mirt() requires named item columns;
ipixuna$items doesn’t have any, so we name them here:
#> item1 item2 item3 item4 item5 item6 item7 item8 item9 item10
#> [1,] 0 1 1 0 1 1 1 1 1 1
#> [2,] 1 0 0 1 0 1 1 1 1 0
#> [3,] 1 1 1 1 0 0 0 1 1 1
#> [4,] 0 0 0 0 0 0 0 0 0 0
#> [5,] 1 0 1 0 1 1 1 1 1 1
#> [6,] 0 0 0 0 0 0 0 0 0 0We fit an exploratory (unrestricted) 2PL item factor model and rotate it (varimax) towards a simple structure:
set.seed(11)
efa <- mirt(items, model = nfactors, itemtype = "2PL", verbose = FALSE)
loadings <- summary(efa, rotate = "varimax", verbose = FALSE)$rotF
loadings#>
#> Loadings:
#> F1 F2 F3
#> item1 0.994
#> item2 0.351 0.437 0.373
#> item3 0.248 0.408 0.655
#> item4 0.741
#> item5 0.967 -0.250
#> item6 0.808 0.150 -0.164
#> item7 0.740 0.115 0.323
#> item8 0.180 0.144 0.750
#> item9 0.588 0.505 0.498
#> item10 0.875 0.240 0.406
#>
#> F1 F2 F3
#> SS loadings 3.469 1.777 2.23Thresholding small loadings to zero turns this into a
0/1 restriction:
#> F1 F2 F3
#> item1 0 1 0
#> item2 1 1 1
#> item3 1 1 1
#> item4 0 0 1
#> item5 1 1 0
#> item6 1 0 0
#> item7 1 0 1
#> item8 0 0 1
#> item9 1 1 1
#> item10 1 0 1This is only one approach, different rotations could be applied
depending on the questionnaire and the latent construct under analysis.
In this example, we will use A to constrain discrimination
parameters to zero.
Step 2: Spatial range with gstat
The second thing vignette("spifa") takes as given is
priors$range, the prior for each Gaussian process’s spatial
range. mirt’s factor scores give a first estimate of each
household’s latent abilities, which we can examine spatially with an
empirical variogram.
Estimate the latent factors
First, we obtain the estimated latent factors using
fscores from mirt:
fscore <- fscores(efa, rotate = "varimax")
colnames(fscore) <- paste0("F", seq_len(nfactors))
head(fscore)#> F1 F2 F3
#> [1,] 1.3211470 -0.08740716 0.40523748
#> [2,] -0.5338226 0.80941335 0.08586527
#> [3,] -0.5509290 1.13178126 1.26640777
#> [4,] -1.1074149 -0.63694759 -0.87035991
#> [5,] 0.9601223 1.01803357 0.05964310
#> [6,] -1.1074149 -0.63694759 -0.87035991Then, we need to geo-reference these factor scores in order to
perform spatial analysis. gstat uses the default
coordinates units to compute distances between locations, so we will
also transform the coordinates system to a planar one where distances
can be computed in meters:
fscore_sf <- st_sf(
id = ipixuna$id, wealth = ipixuna$wealth, as.data.frame(fscore),
geometry = st_transform(st_geometry(ipixuna), 32719)
)
fscore_sf#> Simple feature collection with 100 features and 5 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 201883.2 ymin: 9219014 xmax: 203528.6 ymax: 9221070
#> Projected CRS: WGS 84 / UTM zone 19S
#> First 10 features:
#> id wealth F1 F2 F3 geometry
#> 1 1 -0.57212306 1.32114702 -0.08740716 0.405237478 POINT (202051.7 9219614)
#> 2 2 -0.92003870 -0.53382256 0.80941335 0.085865274 POINT (203142.2 9221070)
#> 3 3 1.23197630 -0.55092905 1.13178126 1.266407773 POINT (203106.4 9220133)
#> 4 4 0.30801774 -1.10741489 -0.63694759 -0.870359908 POINT (202804.6 9220098)
#> 5 5 -0.06247617 0.96012231 1.01803357 0.059643097 POINT (203407.1 9219495)
#> 6 6 2.05414943 -1.10741489 -0.63694759 -0.870359908 POINT (202764.3 9220705)
#> 7 7 2.30552728 -1.42498564 0.88711116 -0.977327784 POINT (202274.5 9219862)
#> 8 8 0.16913811 -0.70727292 0.88660440 -0.009927803 POINT (203219 9220140)
#> 9 9 -0.21989753 0.03349568 -1.44555769 -0.931133875 POINT (202584.1 9219519)
#> 10 10 -0.91744763 0.08459507 1.65413743 0.334133795 POINT (203053.4 9220554)Evaluate predictors effects
We can evaluate the effect of wealth on this estimated latent factors:
fscore_sf |>
tidyr::pivot_longer(F1:F3) |>
ggplot(aes(wealth, value)) +
geom_point() +
geom_smooth(method = 'lm', formula = y ~ x) +
facet_wrap(~ name) +
theme_bw()
There seems to be an effect in F1 and F3, so we will analyse the
residual spatial structure after taking into account wealth
effect:
Evaluate the spatial structure
From gstat, we can use variogram() to
obtain the empirical variogram, vgm() to define a
theoretical variogram and fit.variogram() to fit the
theoretical model to the empirical variogram.
vg1 <- variogram(F1 ~ wealth, fscore_sf, width = 20, cutoff = 1000)
vg1_fit <- fit.variogram(vg1, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg1_fit#> model psill range
#> 1 Nug 0.1913712 0.00000
#> 2 Exp 0.3388614 28.00598
For F1, we can see that there is some evidence of spatial correlation at small scale.
vg2 <- variogram(F2 ~ wealth, fscore_sf, width = 50, cutoff = 1000)
vg2_fit <- fit.variogram(vg2, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg2_fit#> model psill range
#> 1 Nug 0.2078155 0.0000
#> 2 Exp 0.4379639 170.7674
For F2, there is a much clearer spatial structure with a range
parameter around 170.
vg3 <- variogram(F3 ~ wealth, fscore_sf, width = 50, cutoff = 800)
vg3_fit <- fit.variogram(vg3, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg3_fit#> model psill range
#> 1 Nug 0.0000000 0.00000
#> 2 Exp 0.6256558 16.96995
For F3 also the range parameter seems to be small.
In general, there is some evidence of spatial correlation in the
latent factors, and specifically in F2 the spatial structure is strong
with a range parameter around 170. This gives a good understanding of
plausible values for the range parameters of our spifa
model. However, remember that the spifa and
mirt models are not exactly the same, so these values
should be taken only as a reference to decide on our
constraints and priors.
Step 3: Spatial modelling with spifa
With A from Step 1 and a range prior from Step 2 in
hand, spifa() can fit the full spatial model, one
independent Gaussian process per factor by default. For the range prior,
F2’s estimate is the one clear signal – with only 100 households, F1 and
F3’s much smaller ranges are more likely noise than genuine short-range
structure – so we use a common reference value of 150 (close to F2’s
~170) for all three factors:
set.seed(1)
samples <- spifa(
items ~ wealth, data = ipixuna, nfactors = nfactors,
burnin = 1000, niter = 4000,
constraints = list(discrimination = A, loading = diag(nfactors)),
priors = list(range = list(initial = 150, mean = 150, sd = 0.4)))
samples#> Item factor analysis model: spifa_pred
#> Formula: items ~ wealth
#> Dimensions: 100 respondents, 10 items, 3 latent factors, 3 spatial processes
#> MCMC: 1 chain, iter = 4000, thin = 1, samples = 4000
#>
#> Item model parameters:
#> mean median sd q10 q90 ess_bulk rhat
#> c[1] -0.583 -0.554 0.30 -0.974 -0.2272 172 1
#> c[2] -0.581 -0.573 0.25 -0.886 -0.2806 162 1
#> c[3] -0.605 -0.598 0.26 -0.936 -0.2831 166 1
#> c[4] -0.366 -0.353 0.22 -0.654 -0.0958 201 1
#> c[5] -0.330 -0.340 0.27 -0.669 0.0075 208 1
#> c[6] 0.036 0.031 0.22 -0.243 0.3101 264 1
#> c[7] 0.201 0.193 0.24 -0.104 0.5200 148 1
#> c[8] 0.130 0.132 0.22 -0.150 0.4096 216 1
#> c[9] 0.672 0.653 0.34 0.250 1.1239 121 1
#> c[10] 0.398 0.389 0.32 0.011 0.7962 136 1
#> A[2,1] -0.307 -0.300 0.23 -0.599 -0.0170 273 1
#> A[3,1] -0.132 -0.128 0.24 -0.433 0.1608 259 1
#> A[5,1] -1.682 -1.659 0.44 -2.276 -1.1343 114 1
#> A[6,1] -1.228 -1.214 0.32 -1.640 -0.8248 214 1
#> A[7,1] -1.111 -1.082 0.31 -1.514 -0.7368 143 1
#> A[9,1] -1.104 -1.087 0.34 -1.557 -0.6698 124 1
#> A[10,1] -1.401 -1.387 0.36 -1.865 -0.9487 124 1
#> A[1,2] 1.331 1.295 0.43 0.796 1.9100 68 1
#> A[2,2] 0.708 0.690 0.26 0.389 1.0430 183 1
#> A[3,2] 0.599 0.570 0.27 0.281 0.9577 171 1
#> A[5,2] -0.111 -0.105 0.27 -0.457 0.2243 209 1
#> A[9,2] 0.890 0.854 0.38 0.435 1.4033 114 1
#> A[2,3] 0.534 0.528 0.26 0.207 0.8646 275 1
#> A[3,3] 0.946 0.938 0.23 0.653 1.2405 202 1
#> A[4,3] 0.869 0.840 0.29 0.525 1.2438 150 1
#> A[7,3] 0.674 0.662 0.24 0.369 0.9945 285 1
#> A[8,3] 0.972 0.938 0.31 0.617 1.3556 212 1
#> A[9,3] 0.879 0.872 0.30 0.510 1.2577 201 1
#> A[10,3] 1.139 1.107 0.34 0.723 1.6111 151 1
#>
#> Factor model parameters:
#> mean median sd q10 q90 ess_bulk rhat
#> B[1,1] 0.357 0.354 0.106 0.22 0.494 178.6 1.0
#> B[1,2] 0.018 0.013 0.111 -0.12 0.163 327.8 1.0
#> B[1,3] -0.188 -0.185 0.097 -0.31 -0.064 298.4 1.0
#> T[1,1] 0.935 0.925 0.120 0.80 1.106 10.7 1.1
#> T[2,2] 1.014 0.994 0.165 0.83 1.258 26.1 1.0
#> T[3,3] 0.970 0.966 0.148 0.79 1.148 23.6 1.0
#> phi[1] 82.954 80.189 23.631 56.67 112.765 66.2 1.0
#> phi[2] 123.042 120.759 33.020 86.48 169.577 32.1 1.1
#> phi[3] 155.968 154.158 29.957 117.52 198.676 3.3 1.3
#> Corr[2,1] 0.031 0.034 0.422 -0.55 0.630 90.0 1.0
#> Corr[3,1] -0.316 -0.312 0.295 -0.71 0.021 17.2 1.1
#> Corr[3,2] 0.269 0.330 0.420 -0.33 0.782 48.0 1.0
#>
#> ess_bulk is the bulk effective sample size; rhat is the potential
#> scale reduction factor on split chains (Rhat = 1 at convergence).
#> Use summary() for the full set of statistics (incl. ess_tail).Note this uses spifa()’s own default (non-informed)
initial values for everything except the spatial range – the same
starting point a real user, without further domain knowledge, would use.
The full analysis can now proceed using the methods presented in
vignette("spifa") and the package’s other vignettes.
Visualize
Let’s visualize the chains and densities of the posterior. Observe which parameters have better mixing than others.
