Binary CFA

As of version 0.2.1.9000, INLAvaan supports the fitting of binary data for CFA using the pairwise likelihood function available from lavaan. This is an experimental feature which needs further research and testing. Some notes:

• Using PML is considered a “limited-information” approach, due to pairwise likelihood contributions not utilising the fully joint information in the data.
• The scale of the PML function will significantly be different, since the total pairwise log-likelihood adds up contributions in pairs. Following the literature on spatial models using composite likelihoods, INLAvaan adjusts this by a factor of $1/\sqrt{p}$.
• It is proabably unwise to use the Laplace-approximated marginal likelihood for model comparison. The ppp also may not be suitable for ordinal data and needs a rethink.
• Bayesian estimation of CFA favours the parameterization = "theta" option, since the priors on residual variances are more intuitive to specify. INLAvaan switches to this by default.
• For binary models, normal priors for the thresholds should be fine. But looking ahead for ordinal models, the priors for thresholds should be specified in such a way that the ordering is preserved $\tau_0 < \tau_1 < \tau_2 < \cdots < \tau_k$. This is not yet implemented in INLAvaan.

Having said that, let’s take a look at how binary CFA can be implemented in INLAvaan.

library(INLAvaan)
library(blavaan)
#> Loading required package: Rcpp
#> This is blavaan 0.5-10
#> On multicore systems, we suggest use of future::plan("multicore") or
#>   future::plan("multisession") for faster post-MCMC computations.
set.seed(161)

# Generate data
n <- 250
truval <- c(0.8, 0.7, 0.6, 0.5, 0.4, -1.43, -0.55, -0.13, -0.72, -1.13)
dat <- lavaan::simulateData(
  "eta =~ 0.8*y1 + 0.7*5y2 + 0.6*y3 + 0.5*y4 + 0.4*y5
   y1 | -1.43*t1
   y2 | -0.55*t1
   y3 | -0.13*t1
   y4 | -0.72*t1
   y5 | -1.13*t1",
  ordered = TRUE,
  sample.nobs = n
)
head(dat)
#>   y1 y2 y3 y4 y5
#> 1  2  2  1  2  2
#> 2  2  2  2  2  2
#> 3  1  1  2  1  2
#> 4  2  2  2  2  2
#> 5  2  1  1  2  2
#> 6  2  2  2  2  2

# Fit INLAvaan model
mod <- "eta  =~ y1 + y2 + y3 + y4 + y5"
fit <- acfa(mod, dat, ordered = TRUE, std.lv = TRUE, estimator = "PML")
#> ℹ Mode finding and Hessian computation.
#> ✔ Posterior mode and Hessian. [215ms]
#> 
#> ℹ Performing VB correction.
#> ✔ VB correction; mean |δ| = 0.397σ. [522ms]
#> 
#> ⠙ Fitting 0/10 skew-normal marginals.
#> ✔ Fit 10/10 skew-normal marginals. [707ms]
#> 
#> ℹ Adjusting copula correlations (NORTA).
#> ✔ Adjust copula correlations (NORTA). [41ms]
#> 
#> ⠙ Posterior sampling and summarising.
#> ⠹ Computing fit indices (PPP/DIC).
#> ✔ Summarise 1000 posterior draws. [1.2s]
#> 
#> ℹ Fit measures: PPP, DIC.
summary(fit)
#> INLAvaan 0.3.0.9000 ended normally after 37 iterations
#> 
#>   Estimator                                      BAYES
#>   Optimization method                           NLMINB
#>   Number of model parameters                        10
#> 
#>   Number of observations                           250
#> 
#> Model Test (User Model):
#> 
#>    Marginal log-likelihood                   -1125.450 
#>    PPP (Chi-square)                              0.000 
#> 
#> Information Criteria:
#> 
#>    Deviance (DIC)                             2188.439 
#>    Effective parameters (pD)                     7.022 
#> 
#> Parameter Estimates:
#> 
#>    Parameterization                              Theta
#>    Marginalisation method                     SKEWNORM
#>    VB correction                                  TRUE
#> 
#> Latent Variables:
#>                    Estimate       SD     2.5%    97.5%     NMAD    Prior     
#>   eta =~                                                                     
#>     y1                1.109    0.368    0.538    1.954    0.036  normal(0,10)
#>     y2                1.344    0.504    0.611    2.525    0.024  normal(0,10)
#>     y3                0.747    0.249    0.340    1.308    0.008  normal(0,10)
#>     y4                0.272    0.155   -0.011    0.599    0.004  normal(0,10)
#>     y5                0.963    0.321    0.468    1.702    0.019  normal(0,10)
#> 
#> Thresholds:
#>                    Estimate       SD     2.5%    97.5%     NMAD    Prior     
#>     y1|t1            -2.121    0.382   -3.008   -1.547    0.025 normal(0,1.5)
#>     y2|t1            -0.791    0.209   -1.282   -0.491    0.040 normal(0,1.5)
#>     y3|t1            -0.156    0.077   -0.316   -0.012    0.004 normal(0,1.5)
#>     y4|t1            -0.879    0.078   -1.040   -0.734    0.003 normal(0,1.5)
#>     y5|t1            -1.522    0.255   -2.114   -1.141    0.043 normal(0,1.5)
#> 
#> Variances:
#>                    Estimate       SD     2.5%    97.5%     NMAD    Prior     
#>    .y1                1.000                                                  
#>    .y2                1.000                                                  
#>    .y3                1.000                                                  
#>    .y4                1.000                                                  
#>    .y5                1.000                                                  
#>     eta               1.000                                                  
#> 
#> Scales y*:
#>                    Estimate       SD     2.5%    97.5%     NMAD    Prior     
#>     y1                1.519    0.278    1.139    2.154                       
#>     y2                1.694    0.395    1.182    2.626                       
#>     y3                1.261    0.150    1.061    1.629                       
#>     y4                1.046    0.045    1.000    1.172                       
#>     y5                1.400    0.224    1.100    1.952
plot(fit, truth = truval)