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.

Dichotomous Data with Local Dependence and Minor Trait:

  1. Load the package, obtain the data, check the true loading pattern (qlam) and no. of factors/traits.
## LAWBL Package (version 1.4.0; 2021-04-01)
## For tutorials, see https://jinsong-chen.github.io/LAWBL/
dat <- sim24ccfa21$dat
head(dat)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,]    2    1    2    2    1    2    1    1    1     1     1     1     2     1
## [2,]    1    1    1    1    1    1    1    1    1     2     1     1     1     1
## [3,]    1    1    1    1    1    1    1    1    2     2     2     2     1     2
## [4,]    2    2    2    1    2    2    2    2    2     1     1     1     1     1
## [5,]    1    2    1    1    1    1    1    1    1     1     2     1     1     2
## [6,]    1    2    1    1    1    1    2    1    2     1     1     1     2     2
##      [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,]     1     1     2     2     1     2     2     2     1     2
## [2,]     2     1     1     1     1     2     1     2     1     2
## [3,]     2     2     2     2     2     2     2     1     1     2
## [4,]     1     1     2     1     1     2     1     1     2     2
## [5,]     1     1     1     2     1     1     1     1     1     1
## [6,]     2     2     1     1     1     2     1     2     2     2
J <- ncol(dat) # no. of items
qlam <- sim24ccfa21$qlam
qlam
##       [,1] [,2] [,3] [,4]
##  [1,]  0.7  0.0  0.0 0.00
##  [2,]  0.7  0.0  0.0 0.00
##  [3,]  0.7  0.0  0.0 0.00
##  [4,]  0.7  0.0  0.0 0.00
##  [5,]  0.7  0.0  0.0 0.00
##  [6,]  0.7  0.0  0.0 0.00
##  [7,]  0.7  0.0  0.0 0.00
##  [8,]  0.7  0.0  0.0 0.00
##  [9,]  0.0  0.7  0.0 0.00
## [10,]  0.0  0.7  0.0 0.00
## [11,]  0.0  0.7  0.0 0.00
## [12,]  0.0  0.7  0.0 0.00
## [13,]  0.0  0.7  0.0 0.00
## [14,]  0.0  0.7  0.0 0.00
## [15,]  0.0  0.7  0.0 0.55
## [16,]  0.0  0.7  0.0 0.55
## [17,]  0.0  0.0  0.7 0.00
## [18,]  0.0  0.0  0.7 0.00
## [19,]  0.0  0.0  0.7 0.00
## [20,]  0.0  0.0  0.7 0.00
## [21,]  0.0  0.0  0.7 0.00
## [22,]  0.0  0.0  0.7 0.00
## [23,]  0.0  0.0  0.7 0.55
## [24,]  0.0  0.0  0.7 0.55
K <- ncol(qlam) # no. of factors
  1. PCIRM with primary loadings and three traits specified. Longer chain is suggested for stabler performance (burn=iter=5,000 by default).
ipf <- 8
Q<-matrix(-1,J,K-1); # -1 for unspecified items
Q[1:8,1]<-Q[9:16,2]<-Q[17:24,3]<-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
## [19,]   -1   -1    1
## [20,]   -1   -1    1
## [21,]   -1   -1    1
## [22,]   -1   -1    1
## [23,]   -1   -1    1
## [24,]   -1   -1    1
m0<-pcirm(dat = dat,Q=  Q,LD = TRUE,cati = -1,burn = 5000,iter = 5000)
summary(m0)
## $N
## [1] 1000
## 
## $J
## [1] 24
## 
## $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.003180   1.014576
## F2   1.018034   1.022521
## F3   1.004815   1.019459
## 
## $`No. of sig LD terms`
## [1] 7
## 
## $`Cat Items`
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## 
## $`max No. of categories`
## [1] 2
summary(m0, what = 'qlambda')
##          [,1]      [,2]      [,3]
## I1  0.7345649 0.0000000 0.0000000
## I2  0.6973076 0.0000000 0.0000000
## I3  0.6374275 0.0000000 0.0000000
## I4  0.6882978 0.0000000 0.0000000
## I5  0.6594160 0.0000000 0.0000000
## I6  0.7253908 0.0000000 0.0000000
## I7  0.7451140 0.0000000 0.0000000
## I8  0.7600759 0.0000000 0.0000000
## I9  0.0000000 0.7085841 0.0000000
## I10 0.0000000 0.7181708 0.0000000
## I11 0.0000000 0.7004158 0.0000000
## I12 0.0000000 0.6828398 0.0000000
## I13 0.0000000 0.6779277 0.0000000
## I14 0.0000000 0.6817212 0.0000000
## I15 0.0000000 0.7050212 0.0000000
## I16 0.0000000 0.6934798 0.0000000
## I17 0.0000000 0.0000000 0.6843626
## I18 0.0000000 0.0000000 0.7216814
## I19 0.0000000 0.0000000 0.6874567
## I20 0.0000000 0.0000000 0.7218298
## I21 0.0000000 0.0000000 0.7099573
## I22 0.0000000 0.0000000 0.7003339
## I23 0.0000000 0.0000000 0.7146735
## I24 0.0000000 0.0000000 0.7680642
summary(m0, what = 'offpsx') #summarize significant LD terms
##        row col       est         sd      lower     upper sig
## var131   8   7 0.1174882 0.04934009 0.02211229 0.2115044   1
## var247  16  15 0.1581331 0.05443179 0.04929432 0.2610301   1
## var254  23  15 0.1424573 0.04893303 0.04663920 0.2329364   1
## var255  24  15 0.1232736 0.04628765 0.03189747 0.2084133   1
## var263  23  16 0.1268549 0.04819432 0.03587033 0.2190730   1
## var264  24  16 0.1148907 0.04673124 0.02621963 0.2003981   1
## var299  24  23 0.1361735 0.04974018 0.03776283 0.2311301   1
summary(m0,what='int')
##              est         sd       lower      upper sig
## I1   0.002711951 0.02790235 -0.05507659 0.05593188   0
## I2  -0.013572155 0.02850419 -0.06752383 0.04378194   0
## I3   0.002706569 0.03023973 -0.05821248 0.06080019   0
## I4   0.003321130 0.02774261 -0.05065959 0.05697130   0
## I5   0.009994831 0.03005929 -0.04934034 0.06748804   0
## I6  -0.029592488 0.02790604 -0.08449135 0.02543271   0
## I7   0.013743645 0.02759440 -0.04203787 0.06591428   0
## I8   0.015587442 0.02746318 -0.03883204 0.06723259   0
## I9   0.013303222 0.02850328 -0.04013198 0.07144102   0
## I10  0.006658621 0.02871960 -0.04831744 0.06324225   0
## I11 -0.012832268 0.02867042 -0.06903061 0.04285483   0
## I12 -0.013917785 0.02974170 -0.07446867 0.03976769   0
## I13  0.014884230 0.02916929 -0.04095906 0.07365319   0
## I14  0.019664766 0.02905217 -0.03918695 0.07453507   0
## I15 -0.015249240 0.02839058 -0.07069597 0.03989840   0
## I16 -0.022485868 0.02825600 -0.07872775 0.03135618   0
## I17  0.001799232 0.02897047 -0.05621902 0.05736492   0
## I18 -0.019019366 0.02798936 -0.07533539 0.03312753   0
## I19  0.006298121 0.02909733 -0.05234516 0.06063764   0
## I20 -0.015596377 0.02735834 -0.06754596 0.03885849   0
## I21  0.026595115 0.02880197 -0.03054863 0.08137435   0
## I22 -0.013256738 0.02852368 -0.06802863 0.04317556   0
## I23 -0.015661968 0.02758740 -0.06941194 0.03918257   0
## I24 -0.005045578 0.02592744 -0.05493036 0.04558985   0
summary(m0,what='eigen')
##         est        sd    lower    upper sig
## F1 4.077102 0.3004438 3.476225 4.631737   1
## F2 3.976012 0.3075376 3.416388 4.627832   1
## F3 4.174174 0.2918111 3.616153 4.745016   1
#plotting factorial eigenvalue
plot_lawbl(m0) # trace

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

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

  1. PCIRM-LI with four traits and two cross-loadings specified based on results from previous step. Longer chain is suggested for stabler performance.
Q<-cbind(Q,-1);
Q[c(15:16),K] <- 1

m1<-pcirm(dat = dat,Q=  Q,LD = FALSE, cati = -1,burn = 5000,iter = 5000)
summary(m1)
## $N
## [1] 1000
## 
## $J
## [1] 24
## 
## $K
## [1] 4
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] FALSE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 28
## 
## $`True Factor`
## [1] TRUE TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.003113   1.005682
## F2   1.011117   1.018147
## F3   1.007137   1.031432
## F4   1.000765   1.000996
## 
## $`Cat Items`
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## 
## $`max No. of categories`
## [1] 2
summary(m1, what = 'qlambda') #close to qlam
##          [,1]      [,2]      [,3]      [,4]
## I1  0.7343194 0.0000000 0.0000000 0.0000000
## I2  0.7002415 0.0000000 0.0000000 0.0000000
## I3  0.6368656 0.0000000 0.0000000 0.0000000
## I4  0.6845293 0.0000000 0.0000000 0.0000000
## I5  0.6577201 0.0000000 0.0000000 0.0000000
## I6  0.7117158 0.0000000 0.0000000 0.0000000
## I7  0.8236865 0.0000000 0.0000000 0.0000000
## I8  0.8346137 0.0000000 0.0000000 0.0000000
## I9  0.0000000 0.7069743 0.0000000 0.0000000
## I10 0.0000000 0.7382173 0.0000000 0.0000000
## I11 0.0000000 0.7209443 0.0000000 0.0000000
## I12 0.0000000 0.7137747 0.0000000 0.0000000
## I13 0.0000000 0.6908270 0.0000000 0.0000000
## I14 0.0000000 0.6964245 0.0000000 0.0000000
## I15 0.0000000 0.6539522 0.0000000 0.6056571
## I16 0.0000000 0.6424072 0.0000000 0.5672457
## I17 0.0000000 0.0000000 0.6938985 0.0000000
## I18 0.0000000 0.0000000 0.7437660 0.0000000
## I19 0.0000000 0.0000000 0.7038873 0.0000000
## I20 0.0000000 0.0000000 0.7455588 0.0000000
## I21 0.0000000 0.0000000 0.7204857 0.0000000
## I22 0.0000000 0.0000000 0.7077822 0.0000000
## I23 0.0000000 0.0000000 0.6896183 0.4653509
## I24 0.0000000 0.0000000 0.7461177 0.4235851
  1. PCIRM with four traits and four cross-loadings specified based on results from previous step. Longer chain is suggested for stabler performance.
tmp<-summary(m1, what="qlambda")
Q<-matrix(-1,J,K)
Q[tmp!=0]<-1
Q
##       [,1] [,2] [,3] [,4]
##  [1,]    1   -1   -1   -1
##  [2,]    1   -1   -1   -1
##  [3,]    1   -1   -1   -1
##  [4,]    1   -1   -1   -1
##  [5,]    1   -1   -1   -1
##  [6,]    1   -1   -1   -1
##  [7,]    1   -1   -1   -1
##  [8,]    1   -1   -1   -1
##  [9,]   -1    1   -1   -1
## [10,]   -1    1   -1   -1
## [11,]   -1    1   -1   -1
## [12,]   -1    1   -1   -1
## [13,]   -1    1   -1   -1
## [14,]   -1    1   -1   -1
## [15,]   -1    1   -1    1
## [16,]   -1    1   -1    1
## [17,]   -1   -1    1   -1
## [18,]   -1   -1    1   -1
## [19,]   -1   -1    1   -1
## [20,]   -1   -1    1   -1
## [21,]   -1   -1    1   -1
## [22,]   -1   -1    1   -1
## [23,]   -1   -1    1    1
## [24,]   -1   -1    1    1
m2<-pcirm(dat = dat,Q=  Q,LD = TRUE,cati = -1,burn = 5000,iter = 5000)
summary(m2)
## $N
## [1] 1000
## 
## $J
## [1] 24
## 
## $K
## [1] 4
## 
## $`Miss%`
## [1] 0
## 
## $`LD enabled`
## [1] TRUE
## 
## $`Burn in`
## [1] 5000
## 
## $Iteration
## [1] 5000
## 
## $`No. of sig lambda`
## [1] 28
## 
## $`True Factor`
## [1] TRUE TRUE TRUE TRUE
## 
## $`Adj. PSR`
##    Point est. Upper C.I.
## F1   1.008974   1.029809
## F2   1.032254   1.138764
## F3   1.080084   1.313443
## F4   1.110133   1.397390
## 
## $`No. of sig LD terms`
## [1] 1
## 
## $`Cat Items`
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## 
## $`max No. of categories`
## [1] 2
summary(m2, what = 'qlambda')
##          [,1]      [,2]      [,3]      [,4]
## I1  0.7480890 0.0000000 0.0000000 0.0000000
## I2  0.6986163 0.0000000 0.0000000 0.0000000
## I3  0.6410235 0.0000000 0.0000000 0.0000000
## I4  0.6878581 0.0000000 0.0000000 0.0000000
## I5  0.6708255 0.0000000 0.0000000 0.0000000
## I6  0.7136921 0.0000000 0.0000000 0.0000000
## I7  0.7511514 0.0000000 0.0000000 0.0000000
## I8  0.7614955 0.0000000 0.0000000 0.0000000
## I9  0.0000000 0.7018078 0.0000000 0.0000000
## I10 0.0000000 0.7359759 0.0000000 0.0000000
## I11 0.0000000 0.7004909 0.0000000 0.0000000
## I12 0.0000000 0.7017868 0.0000000 0.0000000
## I13 0.0000000 0.6895477 0.0000000 0.0000000
## I14 0.0000000 0.6927472 0.0000000 0.0000000
## I15 0.0000000 0.6349630 0.0000000 0.5569485
## I16 0.0000000 0.6243766 0.0000000 0.5291923
## I17 0.0000000 0.0000000 0.6928552 0.0000000
## I18 0.0000000 0.0000000 0.7228055 0.0000000
## I19 0.0000000 0.0000000 0.6935747 0.0000000
## I20 0.0000000 0.0000000 0.7434658 0.0000000
## I21 0.0000000 0.0000000 0.7020065 0.0000000
## I22 0.0000000 0.0000000 0.6936687 0.0000000
## I23 0.0000000 0.0000000 0.6395868 0.5139557
## I24 0.0000000 0.0000000 0.6943453 0.4722654
summary(m2, what = 'offpsx') #summarize significant LD terms
##        row        col        est         sd      lower      upper        sig 
## 8.00000000 7.00000000 0.11594263 0.04020606 0.03915026 0.19421236 1.00000000
summary(m2,what='int')
##              est         sd       lower      upper sig
## I1   0.002769868 0.02716524 -0.04955226 0.05710893   0
## I2  -0.015144298 0.02852656 -0.07120965 0.03923479   0
## I3   0.002593134 0.02998500 -0.05737031 0.06073508   0
## I4   0.001566732 0.02888058 -0.05759753 0.05621855   0
## I5   0.007415641 0.02941674 -0.05035024 0.06353834   0
## I6  -0.030424656 0.02851169 -0.08441865 0.02689926   0
## I7   0.013617087 0.02729567 -0.03894153 0.06772086   0
## I8   0.014512599 0.02694911 -0.03930752 0.06630339   0
## I9   0.013112777 0.02751001 -0.04102165 0.06641447   0
## I10  0.004911014 0.02771434 -0.04526445 0.06220399   0
## I11 -0.013961139 0.02841734 -0.06730954 0.04352408   0
## I12 -0.015883941 0.02906379 -0.07078718 0.04348562   0
## I13  0.013130438 0.02942875 -0.04691649 0.06899963   0
## I14  0.016523042 0.02810813 -0.03565969 0.07674912   0
## I15 -0.003865385 0.02115651 -0.04309422 0.04006931   0
## I16 -0.010768791 0.02192485 -0.05184711 0.03372420   0
## I17  0.001196895 0.02886823 -0.05576540 0.05641981   0
## I18 -0.020929972 0.02810900 -0.07413061 0.03678940   0
## I19  0.004669439 0.02924214 -0.05454782 0.05943301   0
## I20 -0.017020762 0.02681327 -0.07024862 0.03506574   0
## I21  0.025450933 0.02906586 -0.03169245 0.08174085   0
## I22 -0.013610190 0.02869368 -0.07002762 0.04216997   0
## I23 -0.004888281 0.02115371 -0.04622848 0.03707234   0
## I24  0.003051452 0.02008033 -0.03463350 0.04422939   0
summary(m2,what='eigen')
##         est        sd    lower    upper sig
## F1 4.110620 0.2841569 3.552124 4.632461   1
## F2 3.855238 0.3418707 3.251281 4.550642   1
## F3 3.977654 0.3381574 3.266979 4.619232   1
## F4 1.182787 0.1706071 0.859545 1.531737   1
#plotting factorial eigenvalue
plot_lawbl(m2) # trace

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

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