Note: the estimation process can be time consuming depending on the computing power. You can same some time by reducing the length of the chains.

Continuous Data w/o Local Dependence:

  1. Load the package, obtain the data, and check the true loading pattern (qlam) and local dependence.
## LAWBL Package (version 1.4.0; 2021-04-01)
## For tutorials, see https://jinsong-chen.github.io/LAWBL/
dat <- sim18cfa0$dat
J <- ncol(dat) # no. of items
K <- 3 # no. of factors
sim18cfa0$qlam
##       [,1] [,2] [,3]
##  [1,]  0.7  0.0  0.0
##  [2,]  0.7  0.0  0.0
##  [3,]  0.7  0.0  0.0
##  [4,]  0.7  0.0  0.0
##  [5,]  0.7  0.3  0.0
##  [6,]  0.7  0.3  0.0
##  [7,]  0.0  0.7  0.0
##  [8,]  0.0  0.7  0.0
##  [9,]  0.0  0.7  0.0
## [10,]  0.0  0.7  0.0
## [11,]  0.0  0.7  0.3
## [12,]  0.0  0.7  0.3
## [13,]  0.0  0.0  0.7
## [14,]  0.0  0.0  0.7
## [15,]  0.0  0.0  0.7
## [16,]  0.0  0.0  0.7
## [17,]  0.3  0.0  0.7
## [18,]  0.3  0.0  0.7
sim18cfa0$LD
##      row col
  1. E-step: Estimate with the PCFA-LI model (E-step) by setting LD=F and the design matrix Q. Only a few loadings need to be specified in Q (e.g., 2 per factor). Longer chain is suggested for stabler performance (burn=iter=5,000 by default).
Q<-matrix(-1,J,K); # -1 for unspecified items
Q[1:2,1]<-Q[7:8,2]<-Q[13:14,3]<-1 # 1 for specified items
Q
##       [,1] [,2] [,3]
##  [1,]    1   -1   -1
##  [2,]    1   -1   -1
##  [3,]   -1   -1   -1
##  [4,]   -1   -1   -1
##  [5,]   -1   -1   -1
##  [6,]   -1   -1   -1
##  [7,]   -1    1   -1
##  [8,]   -1    1   -1
##  [9,]   -1   -1   -1
## [10,]   -1   -1   -1
## [11,]   -1   -1   -1
## [12,]   -1   -1   -1
## [13,]   -1   -1    1
## [14,]   -1   -1    1
## [15,]   -1   -1   -1
## [16,]   -1   -1   -1
## [17,]   -1   -1   -1
## [18,]   -1   -1   -1
m0 <- pcfa(dat = dat, Q = Q,LD = FALSE)

# summarize basic information
summary(m0)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] FALSE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.009013   1.013684
## F2   1.020605   1.028176
## F3   1.099059   1.261961
#summarize significant loadings in pattern/Q-matrix format
summary(m0, what = 'qlambda') 
##          [,1]      [,2]      [,3]
## I1  0.7187513 0.0000000 0.0000000
## I2  0.7334514 0.0000000 0.0000000
## I3  0.7395140 0.0000000 0.0000000
## I4  0.7590903 0.0000000 0.0000000
## I5  0.7404959 0.2591514 0.0000000
## I6  0.7131659 0.2757855 0.0000000
## I7  0.0000000 0.7292478 0.0000000
## I8  0.0000000 0.7260071 0.0000000
## I9  0.0000000 0.7375335 0.0000000
## I10 0.0000000 0.7223677 0.0000000
## I11 0.0000000 0.7023648 0.2871526
## I12 0.0000000 0.6923592 0.2685935
## I13 0.0000000 0.0000000 0.7145636
## I14 0.0000000 0.0000000 0.6949948
## I15 0.0000000 0.0000000 0.7275539
## I16 0.0000000 0.0000000 0.7296344
## I17 0.2730959 0.0000000 0.6977168
## I18 0.2605389 0.0000000 0.7021542
#factorial eigenvalue
summary(m0,what='eigen') 
##         est        sd    lower    upper sig
## F1 3.500647 0.5597786 2.634216 4.665791   1
## F2 3.414748 0.4815226 2.533861 4.426002   1
## F3 3.343624 0.5561785 2.260962 4.505174   1
#plotting factorial eigenvalue
plot_lawbl(m0) # trace

plot_lawbl(m0, what='density') #density

plot_lawbl(m0, what='APSR') #adj, PSRF

  1. C-step: Reconfigure the Q matrix for the C-step with one specified loading per item based on results from the E-step. Estimate with the PCFA model by setting LD=TRUE (by default). Longer chain is suggested for stabler performance. Results are very close to the E-step, since there’s no LD in the data.
Q<-matrix(-1,J,K);
tmp<-summary(m0, what="qlambda")
cind<-apply(tmp,1,which.max)
Q[cbind(c(1:J),cind)]<-1
#alternatively
#Q[1:6,1]<-Q[7:12,2]<-Q[13:18,3]<-1 # 1 for specified items

m1 <- pcfa(dat = dat, Q = Q)
summary(m1)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] TRUE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.017644   1.067995
## F2   1.026415   1.098742
## F3   1.009493   1.046284
## 
## $`No. of sig LD terms`
## [1] 0
summary(m1, what = 'qlambda')
##          [,1]      [,2]      [,3]
## I1  0.6774272 0.0000000 0.0000000
## I2  0.7061376 0.0000000 0.0000000
## I3  0.7113447 0.0000000 0.0000000
## I4  0.7300530 0.0000000 0.0000000
## I5  0.7025663 0.2799144 0.0000000
## I6  0.6740022 0.2955424 0.0000000
## I7  0.0000000 0.7354826 0.0000000
## I8  0.0000000 0.7252325 0.0000000
## I9  0.0000000 0.7396601 0.0000000
## I10 0.0000000 0.7244699 0.0000000
## I11 0.0000000 0.7162016 0.2349636
## I12 0.0000000 0.7044436 0.2244517
## I13 0.0000000 0.0000000 0.6967184
## I14 0.0000000 0.0000000 0.6748920
## I15 0.0000000 0.0000000 0.7153157
## I16 0.0000000 0.0000000 0.7051039
## I17 0.2485499 0.0000000 0.6909310
## I18 0.2358830 0.0000000 0.6959562
summary(m1, what = 'offpsx') #summarize significant LD terms
##      row col est sd lower upper sig
summary(m1,what='eigen')
##         est        sd    lower    upper sig
## F1 3.133920 0.2706941 2.616766 3.676457   1
## F2 3.376054 0.3081364 2.818056 3.984023   1
## F3 3.098813 0.2969392 2.546420 3.688709   1
#plotting factorial eigenvalue
# par(mar = rep(2, 4))
plot_lawbl(m1) # trace

plot_lawbl(m1, what='density') #density

plot_lawbl(m1, what='APSR') #adj, PSRF

  1. CFA-LD: One can also configure the Q matrix for a CFA model with local dependence (i.e. without any unspecified loading) based on results from the C-step. Results are also very close.
Q<-summary(m1, what="qlambda")
Q[Q!=0]<-1
Q
##     [,1] [,2] [,3]
## I1     1    0    0
## I2     1    0    0
## I3     1    0    0
## I4     1    0    0
## I5     1    1    0
## I6     1    1    0
## I7     0    1    0
## I8     0    1    0
## I9     0    1    0
## I10    0    1    0
## I11    0    1    1
## I12    0    1    1
## I13    0    0    1
## I14    0    0    1
## I15    0    0    1
## I16    0    0    1
## I17    1    0    1
## I18    1    0    1
m2 <- pcfa(dat = dat, Q = Q)
summary(m2)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] TRUE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.002566   1.009343
## F2   1.006699   1.033092
## F3   1.000344   1.000621
## 
## $`No. of sig LD terms`
## [1] 0
summary(m2, what = 'qlambda') 
##          [,1]      [,2]      [,3]
## I1  0.6831616 0.0000000 0.0000000
## I2  0.6945890 0.0000000 0.0000000
## I3  0.7122678 0.0000000 0.0000000
## I4  0.7097040 0.0000000 0.0000000
## I5  0.6819737 0.3297539 0.0000000
## I6  0.6690497 0.3418172 0.0000000
## I7  0.0000000 0.7293637 0.0000000
## I8  0.0000000 0.7234732 0.0000000
## I9  0.0000000 0.7232783 0.0000000
## I10 0.0000000 0.6929052 0.0000000
## I11 0.0000000 0.6918715 0.2902011
## I12 0.0000000 0.6899762 0.2915237
## I13 0.0000000 0.0000000 0.6759039
## I14 0.0000000 0.0000000 0.6797780
## I15 0.0000000 0.0000000 0.7212977
## I16 0.0000000 0.0000000 0.6672593
## I17 0.3107708 0.0000000 0.6791115
## I18 0.2880834 0.0000000 0.6816405
summary(m2, what = 'offpsx')
##      row col est sd lower upper sig
summary(m2,what='eigen')
##         est        sd    lower    upper sig
## F1 3.062037 0.1659698 2.743449 3.393940   1
## F2 3.247821 0.1784674 2.905209 3.589116   1
## F3 2.989231 0.1671173 2.664989 3.315374   1
plot_lawbl(m2) # Eigens' traces are excellent without regularization of the loadings

Continuous Data with Local Dependence:

  1. Obtain the data and check the true loading pattern (qlam) and local dependence.
dat <- sim18cfa1$dat
J <- ncol(dat) # no. of items
K <- 3 # no. of factors
sim18cfa1$qlam
##       [,1] [,2] [,3]
##  [1,]  0.7  0.0  0.0
##  [2,]  0.7  0.0  0.0
##  [3,]  0.7  0.0  0.0
##  [4,]  0.7  0.0  0.0
##  [5,]  0.7  0.3  0.0
##  [6,]  0.7  0.3  0.0
##  [7,]  0.0  0.7  0.0
##  [8,]  0.0  0.7  0.0
##  [9,]  0.0  0.7  0.0
## [10,]  0.0  0.7  0.0
## [11,]  0.0  0.7  0.3
## [12,]  0.0  0.7  0.3
## [13,]  0.0  0.0  0.7
## [14,]  0.0  0.0  0.7
## [15,]  0.0  0.0  0.7
## [16,]  0.0  0.0  0.7
## [17,]  0.3  0.0  0.7
## [18,]  0.3  0.0  0.7
sim18cfa1$LD # effect size = .3
##      row col
## [1,]  14   1
## [2,]   7   2
## [3,]   4   3
## [4,]  13   8
## [5,]  10   9
## [6,]  16  15
  1. E-step: Estimate with the PCFA-LI model (E-step) by setting LD=FALSE and the design matrix Q. Only a few loadings need to be specified in Q (e.g., 2 per factor). Some loading estimates are biased due to ignoring the LD. So do the eigenvalues.
Q<-matrix(-1,J,K); # -1 for unspecified items
Q[1:2,1]<-Q[7:8,2]<-Q[13:14,3]<-1 # 1 for specified items
Q
##       [,1] [,2] [,3]
##  [1,]    1   -1   -1
##  [2,]    1   -1   -1
##  [3,]   -1   -1   -1
##  [4,]   -1   -1   -1
##  [5,]   -1   -1   -1
##  [6,]   -1   -1   -1
##  [7,]   -1    1   -1
##  [8,]   -1    1   -1
##  [9,]   -1   -1   -1
## [10,]   -1   -1   -1
## [11,]   -1   -1   -1
## [12,]   -1   -1   -1
## [13,]   -1   -1    1
## [14,]   -1   -1    1
## [15,]   -1   -1   -1
## [16,]   -1   -1   -1
## [17,]   -1   -1   -1
## [18,]   -1   -1   -1
m0 <- pcfa(dat = dat, Q = Q,LD = FALSE)
summary(m0)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] FALSE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.047566   1.108558
## F2   1.397283   2.373383
## F3   1.039068   1.087543
summary(m0, what = 'qlambda')
##          [,1]      [,2]      [,3]
## I1  0.6218244 0.0000000 0.0000000
## I2  0.6760114 0.0000000 0.0000000
## I3  0.8672665 0.0000000 0.0000000
## I4  0.8568885 0.0000000 0.0000000
## I5  0.6937996 0.2609903 0.0000000
## I6  0.6974870 0.2645690 0.0000000
## I7  0.0000000 0.6319171 0.0000000
## I8  0.0000000 0.6542019 0.0000000
## I9  0.0000000 0.8667263 0.0000000
## I10 0.0000000 0.8838674 0.0000000
## I11 0.0000000 0.6350423 0.2956779
## I12 0.0000000 0.6415892 0.2723378
## I13 0.0000000 0.0000000 0.7267767
## I14 0.0000000 0.0000000 0.6913619
## I15 0.0000000 0.0000000 0.7946597
## I16 0.0000000 0.0000000 0.8069249
## I17 0.2952332 0.0000000 0.6829470
## I18 0.2954748 0.0000000 0.6876156
summary(m0,what='eigen')
##         est        sd    lower    upper sig
## F1 3.598526 0.5254012 2.722497 4.555573   1
## F2 3.452701 0.5762372 2.497630 4.541254   1
## F3 3.583772 0.4327281 2.788800 4.417996   1
plot_lawbl(m0) # trace

plot_lawbl(m0, what='APSR')

  1. C-step: Reconfigure the Q matrix for the C-step with one specified loading per item based on results from the E-step. Estimate with the PCFA model by setting LD=TRUE (by default). The estimates are more accurate, and the LD terms can be largely recovered.
Q<-matrix(-1,J,K);
tmp<-summary(m0, what="qlambda")
cind<-apply(tmp,1,which.max)
Q[cbind(c(1:J),cind)]<-1
Q
##       [,1] [,2] [,3]
##  [1,]    1   -1   -1
##  [2,]    1   -1   -1
##  [3,]    1   -1   -1
##  [4,]    1   -1   -1
##  [5,]    1   -1   -1
##  [6,]    1   -1   -1
##  [7,]   -1    1   -1
##  [8,]   -1    1   -1
##  [9,]   -1    1   -1
## [10,]   -1    1   -1
## [11,]   -1    1   -1
## [12,]   -1    1   -1
## [13,]   -1   -1    1
## [14,]   -1   -1    1
## [15,]   -1   -1    1
## [16,]   -1   -1    1
## [17,]   -1   -1    1
## [18,]   -1   -1    1
m1 <- pcfa(dat = dat, Q = Q)
summary(m1)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] TRUE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.046817   1.105924
## F2   1.109160   1.403220
## F3   1.032252   1.103626
## 
## $`No. of sig LD terms`
## [1] 6
summary(m1, what = 'qlambda')
##          [,1]      [,2]      [,3]
## I1  0.6768573 0.0000000 0.0000000
## I2  0.6889138 0.0000000 0.0000000
## I3  0.7245866 0.0000000 0.0000000
## I4  0.7132980 0.0000000 0.0000000
## I5  0.7422758 0.2227292 0.0000000
## I6  0.7262931 0.2271932 0.0000000
## I7  0.0000000 0.7049859 0.0000000
## I8  0.0000000 0.6932314 0.0000000
## I9  0.0000000 0.7147083 0.0000000
## I10 0.0000000 0.7298097 0.0000000
## I11 0.0000000 0.7069986 0.2413025
## I12 0.0000000 0.7176572 0.2165244
## I13 0.0000000 0.0000000 0.7658746
## I14 0.0000000 0.0000000 0.6966367
## I15 0.0000000 0.0000000 0.7201338
## I16 0.0000000 0.0000000 0.7327813
## I17 0.2329523 0.0000000 0.7168271
## I18 0.2344990 0.0000000 0.7222813
summary(m1,what='eigen')
##         est        sd    lower    upper sig
## F1 3.252581 0.3281683 2.670066 3.929519   1
## F2 3.242464 0.3922696 2.568701 4.086261   1
## F3 3.374163 0.3300497 2.759612 4.049376   1
summary(m1, what = 'offpsx')
##        row col       est         sd     lower     upper sig
## var14   14   1 0.2636372 0.04955688 0.1571931 0.3548216   1
## var24    7   2 0.2505906 0.04790738 0.1483714 0.3428096   1
## var37    4   3 0.3037038 0.06701136 0.1799548 0.4420304   1
## var111  13   8 0.2328172 0.04150166 0.1420081 0.3095174   1
## var118  10   9 0.3064860 0.06224591 0.1807960 0.4254826   1
## var163  16  15 0.2730498 0.06066370 0.1437833 0.3925561   1
  1. CFA-LD: Configure the Q matrix for a CFA model with local dependence (i.e. without any unspecified loading) based on results from the C-step. Results are better than, but similar to the C-step.
Q<-summary(m1, what="qlambda")
Q[Q!=0]<-1
Q
##     [,1] [,2] [,3]
## I1     1    0    0
## I2     1    0    0
## I3     1    0    0
## I4     1    0    0
## I5     1    1    0
## I6     1    1    0
## I7     0    1    0
## I8     0    1    0
## I9     0    1    0
## I10    0    1    0
## I11    0    1    1
## I12    0    1    1
## I13    0    0    1
## I14    0    0    1
## I15    0    0    1
## I16    0    0    1
## I17    1    0    1
## I18    1    0    1
m2 <- pcfa(dat = dat, Q = Q)
summary(m2)
## $N
## [1] 1000
## 
## $J
## [1] 18
## 
## $K
## [1] 3
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] TRUE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 24
## 
## $`True Factor`
## [1] TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.001739   1.005549
## F2   1.025481   1.112992
## F3   1.007960   1.024105
## 
## $`No. of sig LD terms`
## [1] 6
summary(m2, what = 'qlambda') 
##          [,1]      [,2]      [,3]
## I1  0.6926975 0.0000000 0.0000000
## I2  0.6749251 0.0000000 0.0000000
## I3  0.7290900 0.0000000 0.0000000
## I4  0.7146251 0.0000000 0.0000000
## I5  0.7103587 0.3031598 0.0000000
## I6  0.6978269 0.3117567 0.0000000
## I7  0.0000000 0.6922642 0.0000000
## I8  0.0000000 0.7120908 0.0000000
## I9  0.0000000 0.6836637 0.0000000
## I10 0.0000000 0.6805879 0.0000000
## I11 0.0000000 0.6821991 0.3235585
## I12 0.0000000 0.6889061 0.2949464
## I13 0.0000000 0.0000000 0.7155813
## I14 0.0000000 0.0000000 0.6828643
## I15 0.0000000 0.0000000 0.7286837
## I16 0.0000000 0.0000000 0.7363161
## I17 0.3043890 0.0000000 0.7083721
## I18 0.2995753 0.0000000 0.6916461
summary(m2,what='eigen')
##         est        sd    lower    upper sig
## F1 3.162347 0.1762245 2.837415 3.512905   1
## F2 3.057030 0.1722605 2.732651 3.402668   1
## F3 3.234205 0.1890688 2.878885 3.614345   1
summary(m2, what = 'offpsx')
##        row col       est         sd     lower     upper sig
## var14   14   1 0.2914881 0.02921900 0.2347119 0.3477505   1
## var24    7   2 0.2900500 0.03020686 0.2304911 0.3502030   1
## var37    4   3 0.2716046 0.04984652 0.1738372 0.3709615   1
## var111  13   8 0.2752596 0.03022265 0.2192931 0.3354367   1
## var118  10   9 0.3246022 0.04956960 0.2323996 0.4246607   1
## var163  16  15 0.2561622 0.04685895 0.1612006 0.3394688   1