Advanced features

This vignette goes deeper into two things vignette("spifa") and vignette("vg02-model-types") only touch briefly:

  • how to control constraints/priors for every parameter group,
  • how to continue sampling with update() when a fit hasn’t converged yet, and
  • how to compare models using dic().
library(spifa)

data(ipixuna)
nfactors <- 3
nitems <- ncol(ipixuna$items)

Restrictions and priors

spifa()’s parameter glossary (see ?spifa) maps six configurable parameter names – easiness (c), discrimination (A), effect (B), corr (Corr), loading (T), range (phi). Two spifa() arguments use these names:

  • constraints: fixes cells of a parameter to a value instead of estimating them. Only discrimination, loading, and sd support this (the model needs some fixed structure for identifiability).
  • priors: sets the initial value and prior hyperparameters (mean/sd, or eta for corr) for a parameter that is estimated. Every one of the six names accepts this.

constraints

  • discrimination (nitems x nfactors matrix): which items inform which factors,
  • loading (nfactors x ngp matrix): which Gaussian processes feed into which factors, and
  • sd: fixed values for the latent factors’ residual standard deviation (length nfactors). Unlike discrimination/loading, this is not estimated, it is always fixed.

Define discrimination structure:

A <- matrix(1, nitems, nfactors)
A[c(4, 8), 1] <- 0
A[c(4, 5, 6, 7, 8, 10), 2] <- 0
A[c(5, 6), 3] <- 0
A
#>       [,1] [,2] [,3]
#>  [1,]    1    1    1
#>  [2,]    1    1    1
#>  [3,]    1    1    1
#>  [4,]    0    0    1
#>  [5,]    1    0    0
#>  [6,]    1    0    0
#>  [7,]    1    0    1
#>  [8,]    0    0    1
#>  [9,]    1    1    1
#> [10,]    1    0    1

Define loading structure, factors 1 and 2 share one Gaussian process, factor 3 gets its own:

ngp <- 2
T_loading <- matrix(c(1, 1, 0, 0, 0, 1), nfactors, ngp)
T_loading
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    1    0
#> [3,]    0    1

Fix the residual standard deviation of each factor to 0.2:

sd_fixed <- rep(0.2, nfactors)
sd_fixed
#> [1] 0.2 0.2 0.2

priors

Every element accepts initial (starting value for the sampler) and the prior hyperparameters mean/sd. The parameter corr is the exception: it takes initial and eta (the LKJ shape parameter) instead.

By default, spifa() starts every parameter at a generic default and lets the data speak. If you have domain knowledge, setting priors closer to a plausible value can speed up convergence substantially. Here we set a plausible prior for every parameter group:

  • item easiness roughly increasing across items,
  • a couple of items with a strong sign-fixing discrimination prior (as in vignette("spifa")),
  • a wealth effect that differs in sign across factors,
  • moderate loadings,
  • factor correlations nudged away from independence, and
  • a range parameter around 150 as in vignette("spifa"):
# easiness hyperparameters
c_mean <- seq(-0.7, 0.9, length.out = nitems)
c_sd <- rep(0.3, nitems)

# discrimination hyperparameters
A_mean <- matrix(0, nitems, nfactors)
A_mean[c(5, 6), 1] <- 1
A_mean[c(4, 8), 3] <- -1
A_sd <- matrix(1, nitems, nfactors)
A_sd[A_mean != 0] <- 0.45

# predictor effects hyperparameters
B_mean <- matrix(c(-0.5, 0.5, 0.2), nrow = 1)
B_sd <- matrix(0.5, nrow = 1, ncol = nfactors)

# correlation initial value
corr_initial <- diag(nfactors)

# range parameter hyperparameters
phi_mean <- 150
phi_sd <- 0.4

Fit with constraints and priors

In the following code, we explicitely defined the constraints, and priors for all possible parameters:

set.seed(1)
fit_informed <- spifa(
  items ~ wealth, data = ipixuna, nfactors = nfactors, ngp = 3,
  burnin = 200, niter = 200, thin = 2,
  constraints = list(discrimination = A, loading = diag(nfactors), sd = sd_fixed),
  priors = list(
    easiness = list(initial = c_mean, mean = c_mean, sd = c_sd),
    discrimination = list(initial = A_mean, mean = A_mean, sd = A_sd),
    effect = list(initial = B_mean, mean = B_mean, sd = B_sd),
    corr = list(initial = corr_initial, eta = 2),
    loading = list(initial = 0.6, mean = 0.6, sd = 0.3),
    range = list(initial = phi_mean, mean = phi_mean, sd = phi_sd)))
fit_informed
#> Item factor analysis model: spifa_pred
#> Formula: items ~ wealth
#> Dimensions: 100 respondents, 10 items, 3 latent factors, 3 spatial processes
#> MCMC: 1 chain, iter = 199, thin = 2, samples = 100
#> 
#> Item model parameters:
#>           mean median   sd    q10    q90 ess_bulk rhat
#> c[1]    -0.336 -0.333 0.19 -0.553 -0.088     24.5  1.1
#> c[2]    -0.358 -0.368 0.18 -0.563 -0.119      4.2  1.3
#> c[3]    -0.445 -0.459 0.21 -0.707 -0.197     11.6  1.1
#> c[4]    -0.285 -0.278 0.14 -0.473 -0.101     29.7  1.0
#> c[5]    -0.283 -0.296 0.17 -0.490 -0.065     44.0  1.0
#> c[6]    -0.038 -0.053 0.13 -0.188  0.116      7.3  1.2
#> c[7]     0.093  0.078 0.16 -0.114  0.302      3.9  1.2
#> c[8]     0.100  0.083 0.15 -0.076  0.319      5.8  1.2
#> c[9]     0.486  0.488 0.19  0.234  0.736     23.1  1.0
#> c[10]    0.223  0.182 0.22 -0.028  0.508      7.2  1.1
#> A[1,1]   0.724  0.696 0.31  0.362  1.065     12.8  1.0
#> A[2,1]   0.408  0.359 0.24  0.136  0.735      3.0  1.3
#> A[3,1]   0.277  0.273 0.17  0.072  0.506     30.8  1.1
#> A[5,1]   1.239  1.231 0.24  0.934  1.547     32.1  1.0
#> A[6,1]   1.059  1.055 0.20  0.807  1.298     21.7  1.1
#> A[7,1]   1.119  1.124 0.30  0.729  1.504      2.9  1.3
#> A[9,1]   1.212  1.120 0.39  0.759  1.808      2.3  1.4
#> A[10,1]  1.349  1.318 0.24  1.073  1.667     28.4  1.0
#> A[1,2]   1.368  1.320 0.33  0.984  1.807      2.6  1.4
#> A[2,2]   0.496  0.489 0.23  0.230  0.778      7.0  1.1
#> A[3,2]   0.412  0.392 0.24  0.108  0.735     27.4  1.0
#> A[9,2]   0.608  0.575 0.30  0.266  1.024      2.2  1.4
#> A[1,3]  -0.463 -0.466 0.22 -0.732 -0.222     14.3  1.0
#> A[2,3]  -0.701 -0.657 0.22 -0.988 -0.485     29.7  1.0
#> A[3,3]  -1.206 -1.182 0.32 -1.626 -0.863     10.4  1.1
#> A[4,3]  -0.774 -0.770 0.21 -1.047 -0.506     23.6  1.0
#> A[7,3]  -0.500 -0.484 0.17 -0.705 -0.280     22.5  1.0
#> A[8,3]  -0.871 -0.877 0.20 -1.125 -0.608     25.4  1.1
#> A[9,3]  -0.853 -0.839 0.24 -1.193 -0.520      2.0  1.5
#> A[10,3] -1.015 -1.019 0.27 -1.354 -0.634     20.4  1.0
#> 
#> Factor model parameters:
#>              mean  median     sd     q10    q90 ess_bulk rhat
#> B[1,1]     -0.682  -0.678  0.077  -0.788  -0.58     69.6 1.01
#> B[1,2]      0.671   0.669  0.077   0.581   0.78     99.6 0.99
#> B[1,3]      0.342   0.350  0.083   0.242   0.44     75.0 1.00
#> T[1,1]      0.728   0.748  0.033   0.668   0.76      7.4 1.11
#> T[2,2]      0.719   0.718  0.041   0.649   0.75     13.0 1.03
#> T[3,3]      0.923   0.922  0.069   0.827   1.01      4.1 1.29
#> phi[1]    148.628 145.465 21.836 127.235 172.44     12.1 1.11
#> phi[2]    146.275 138.778 22.228 120.307 175.42      5.0 1.22
#> phi[3]    125.498 118.301 21.699 100.789 160.64     11.6 1.02
#> Corr[2,1]   0.072   0.083  0.072   0.010   0.14      5.6 1.16
#> Corr[3,1]   0.100   0.105  0.092  -0.063   0.19     13.6 1.04
#> Corr[3,2]   0.098   0.113  0.094  -0.078   0.20     10.1 1.09
#> 
#> 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).

Continue sampling

Fits don’t always converge on the first try. If you take a look at rhat in the previous section, you will notice that some values are much higher that 1.05 for example. To continue sampling from the current state you can use update():

fit1 <- update(fit_informed, niter = 1000)
fit1
#> Item factor analysis model: spifa_pred
#> Formula: items ~ wealth
#> Dimensions: 100 respondents, 10 items, 3 latent factors, 3 spatial processes
#> MCMC: 1 chain, iter = 1000, thin = 1, samples = 1000
#> 
#> Item model parameters:
#>            mean median   sd    q10    q90 ess_bulk rhat
#> c[1]    -0.2638 -0.260 0.18 -0.493 -0.043    135.5  1.0
#> c[2]    -0.3776 -0.363 0.18 -0.620 -0.160     76.1  1.0
#> c[3]    -0.4079 -0.409 0.19 -0.635 -0.159    100.3  1.0
#> c[4]    -0.2783 -0.287 0.16 -0.483 -0.062     61.4  1.0
#> c[5]    -0.2751 -0.279 0.17 -0.497 -0.048    128.8  1.0
#> c[6]     0.0031  0.011 0.18 -0.222  0.234    183.8  1.0
#> c[7]     0.1404  0.142 0.17 -0.074  0.349    135.0  1.0
#> c[8]     0.1327  0.124 0.17 -0.074  0.367    122.2  1.0
#> c[9]     0.4895  0.500 0.19  0.237  0.712    107.1  1.0
#> c[10]    0.2780  0.275 0.21  0.019  0.534     82.7  1.0
#> A[1,1]   0.5636  0.546 0.32  0.182  0.993     21.7  1.1
#> A[2,1]   0.5803  0.564 0.30  0.203  0.959     39.3  1.0
#> A[3,1]   0.2813  0.294 0.27 -0.073  0.612     72.4  1.0
#> A[5,1]   1.3450  1.326 0.26  1.031  1.677     58.0  1.0
#> A[6,1]   1.3215  1.311 0.26  0.995  1.661     43.0  1.1
#> A[7,1]   1.0145  1.011 0.26  0.691  1.342     65.4  1.0
#> A[9,1]   1.1041  1.079 0.37  0.645  1.589     43.4  1.0
#> A[10,1]  1.4137  1.374 0.40  0.933  1.961     28.2  1.0
#> A[1,2]   1.1244  1.086 0.38  0.645  1.638     36.5  1.0
#> A[2,2]   0.7975  0.751 0.34  0.409  1.269     34.7  1.0
#> A[3,2]   0.4960  0.493 0.29  0.152  0.863     60.4  1.0
#> A[9,2]   0.5556  0.520 0.31  0.185  0.978     63.7  1.0
#> A[1,3]  -0.5143 -0.497 0.31 -0.914 -0.153     51.2  1.0
#> A[2,3]  -0.8180 -0.796 0.27 -1.164 -0.488     41.9  1.0
#> A[3,3]  -1.2694 -1.236 0.34 -1.703 -0.862     38.8  1.1
#> A[4,3]  -0.8805 -0.867 0.27 -1.255 -0.528     33.7  1.0
#> A[7,3]  -0.6520 -0.634 0.24 -0.959 -0.362     75.6  1.0
#> A[8,3]  -1.0112 -1.004 0.25 -1.308 -0.715     47.9  1.0
#> A[9,3]  -1.0773 -1.040 0.35 -1.522 -0.665     53.9  1.0
#> A[10,3] -1.2777 -1.232 0.44 -1.861 -0.769      4.9  1.2
#> 
#> Factor model parameters:
#>              mean  median     sd    q10    q90 ess_bulk rhat
#> B[1,1]     -0.637  -0.636  0.087  -0.75  -0.53      374  1.0
#> B[1,2]      0.730   0.727  0.070   0.64   0.82      672  1.0
#> B[1,3]      0.437   0.432  0.095   0.32   0.56      193  1.0
#> T[1,1]      0.776   0.787  0.062   0.68   0.85       33  1.0
#> T[2,2]      0.639   0.636  0.048   0.57   0.71       28  1.1
#> T[3,3]      0.847   0.847  0.140   0.66   1.03       13  1.0
#> phi[1]    109.023 104.894 33.715  67.87 153.61       16  1.0
#> phi[2]    190.876 185.688 32.882 159.08 241.64       24  1.1
#> phi[3]    146.985 141.984 40.638  99.85 217.74       32  1.0
#> Corr[2,1]   0.246   0.233  0.080   0.16   0.36       45  1.0
#> Corr[3,1]   0.020   0.026  0.178  -0.22   0.25       19  1.1
#> Corr[3,2]  -0.027  -0.017  0.194  -0.26   0.22       12  1.1
#> 
#> 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).

rhat has moved closer to 1 with the extra iterations. In practice, you have to run many more iterations.

Compare models with DIC

dic() computes the Deviance Information Criterion for a fitted model, useful for comparing candidate models fitted to the same data. You should compare them once convergence is ensured. For example, let’s compare with a spifa with 2 latent factors and default constraints and priors:

fit2 <- spifa(items ~ wealth, data = ipixuna, nfactors = 2, niter = 1000, burnin = 1000)
fit2
#> Item factor analysis model: spifa_pred
#> Formula: items ~ wealth
#> Dimensions: 100 respondents, 10 items, 2 latent factors, 2 spatial processes
#> MCMC: 1 chain, iter = 1000, thin = 1, samples = 1000
#> 
#> Item model parameters:
#>           mean median   sd     q10    q90 ess_bulk rhat
#> c[1]    -0.383 -0.380 0.14 -0.5612 -0.207    185.5  1.0
#> c[2]    -0.507 -0.508 0.17 -0.7356 -0.293    145.3  1.0
#> c[3]    -0.628 -0.592 0.26 -0.9625 -0.330     34.5  1.0
#> c[4]    -0.324 -0.325 0.17 -0.5443 -0.112    183.6  1.0
#> c[5]    -0.332 -0.320 0.24 -0.6406 -0.035     55.2  1.0
#> c[6]     0.051  0.044 0.19 -0.1995  0.303    152.2  1.0
#> c[7]     0.226  0.225 0.18 -0.0011  0.467    124.4  1.0
#> c[8]     0.193  0.187 0.19 -0.0438  0.446    105.6  1.0
#> c[9]     0.722  0.704 0.29  0.3456  1.110     37.0  1.0
#> c[10]    0.500  0.492 0.28  0.1602  0.849     72.0  1.0
#> A[1,1]   0.420  0.408 0.18  0.2019  0.684     90.6  1.0
#> A[2,1]   0.655  0.642 0.29  0.3019  1.014     24.8  1.0
#> A[3,1]   0.817  0.798 0.45  0.2214  1.424     10.2  1.1
#> A[4,1]   0.296  0.238 0.41 -0.1768  0.882     11.8  1.0
#> A[5,1]   1.666  1.610 0.60  0.9255  2.437      7.7  1.1
#> A[6,1]   1.254  1.250 0.41  0.7358  1.784     19.3  1.0
#> A[7,1]   1.394  1.359 0.40  0.9059  1.950     45.6  1.0
#> A[8,1]   0.568  0.565 0.37  0.0861  1.031     12.2  1.1
#> A[9,1]   1.691  1.670 0.48  1.0592  2.312     36.3  1.0
#> A[10,1]  2.065  2.042 0.55  1.4238  2.713     17.6  1.0
#> A[2,2]   0.538  0.519 0.26  0.2131  0.891     28.8  1.0
#> A[3,2]   0.987  0.915 0.48  0.4372  1.661     14.5  1.1
#> A[4,2]   0.891  0.830 0.38  0.4779  1.389     13.8  1.1
#> A[5,2]  -0.910 -0.936 0.51 -1.5638 -0.226     12.5  1.0
#> A[6,2]  -0.910 -0.883 0.49 -1.5376 -0.315      9.2  1.1
#> A[7,2]  -0.097 -0.071 0.43 -0.7189  0.438     10.1  1.1
#> A[8,2]   0.846  0.825 0.35  0.3936  1.284     29.2  1.0
#> A[9,2]   0.435  0.363 0.54 -0.1897  1.176      9.5  1.1
#> A[10,2]  0.200  0.178 0.54 -0.4857  0.904      8.9  1.1
#> 
#> Factor model parameters:
#>              mean  median     sd    q10    q90 ess_bulk rhat
#> B[1,1]    -0.4397 -0.4392 0.1172 -0.584 -0.297      341  1.0
#> B[1,2]    -0.0018 -0.0073 0.2069 -0.273  0.284       11  1.1
#> T[1,1]     0.4976  0.4795 0.0871  0.399  0.662       13  1.1
#> T[2,2]     0.5316  0.5252 0.0834  0.420  0.627       30  1.0
#> phi[1]     0.0488  0.0484 0.0100  0.036  0.064       27  1.0
#> phi[2]     0.0515  0.0504 0.0094  0.040  0.066       52  1.0
#> Corr[2,1]  0.2155  0.2392 0.4061 -0.437  0.738        9  1.2
#> 
#> 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).

Computing dic() for each model:

# M1: 3 factors
dic(fit1, burnin = 100, thin = 2)
#> # A tibble: 1 × 3
#>   mean_deviance p_eff   dic
#>           <dbl> <dbl> <dbl>
#> 1          850.  148.  998.
# M2: 2 factors
dic(fit2, burnin = 100, thin = 2)
#> # A tibble: 1 × 3
#>   mean_deviance p_eff   dic
#>           <dbl> <dbl> <dbl>
#> 1          880.  147. 1027.

Assuming both models have converged, the 3-factor model has the lower DIC, favouring the discrimination structure used above over an 2-factor alternative.