library(blcfa)
## Loading required package: stringr
## Loading required package: MASS
## Loading required package: MCMCpack
## Loading required package: coda
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2024 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
## Loading required package: msm
## Loading required package: statmod
## Loading required package: psychometric
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: multilevel
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: purrr
## Loading required package: MplusAutomation
## Version: 1.1.1
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
## Loading required package: sna
## Loading required package: statnet.common
##
## Attaching package: 'statnet.common'
## The following objects are masked from 'package:base':
##
## attr, order
## Loading required package: network
##
## 'network' 1.18.2 (2023-12-04), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information
## sna: Tools for Social Network Analysis
## Version 2.7-2 created on 2023-12-05.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## For citation information, type citation("sna").
## Type help(package="sna") to get started.
##
## Attaching package: 'sna'
## The following object is masked from 'package:nlme':
##
## gapply
## Loading required package: doParallel
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: rstudioapi
library(blavaan)
## Loading required package: Rcpp
## This is blavaan 0.5-4
## On multicore systems, we suggest use of future::plan("multicore") or
## future::plan("multisession") for faster post-MCMC computations.
library(foreign)
library(plyr)
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following object is masked from 'package:network':
##
## is.discrete
## The following object is masked from 'package:purrr':
##
## compact
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
library(blavaan)
library(semTools)
## Loading required package: lavaan
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
##
## ###############################################################################
## This is semTools 0.5-6
## All users of R (or SEM) are invited to submit functions or ideas for functions.
## ###############################################################################
library(rstan)
## Loading required package: StanHeaders
##
## rstan version 2.32.6 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
##
## Attaching package: 'rstan'
## The following object is masked from 'package:MplusAutomation':
##
## extract
## The following object is masked from 'package:coda':
##
## traceplot
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:psychometric':
##
## alpha
library(HDInterval)
library(mcmcse)
library(Rmpfr)
## Loading required package: gmp
##
## Attaching package: 'gmp'
## The following objects are masked from 'package:base':
##
## %*%, apply, crossprod, matrix, tcrossprod
## C code of R package 'Rmpfr': GMP using 64 bits per limb
##
## Attaching package: 'Rmpfr'
## The following object is masked from 'package:gmp':
##
## outer
## The following objects are masked from 'package:stats':
##
## dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm
## The following objects are masked from 'package:base':
##
## cbind, pmax, pmin, rbind
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:psychometric':
##
## alpha
## The following object is masked from 'package:purrr':
##
## discard
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(epiDisplay)
## Loading required package: survival
## Loading required package: nnet
##
## Attaching package: 'epiDisplay'
## The following object is masked from 'package:scales':
##
## alpha
## The following object is masked from 'package:ggplot2':
##
## alpha
## The following object is masked from 'package:rstan':
##
## lookup
## The following object is masked from 'package:psychometric':
##
## alpha
library(bayesplot)
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
library(Bayesrel)
RIASEC42 = read.csv2("D:\\OneDrive - UMP\\R - lenh 2016\\RIASEC42.csv", comment.char="#", stringsAsFactors=T)
names(RIASEC42)
## [1] "id" "R1" "R2" "R3" "R4" "R5" "R6" "R7" "I1" "I2" "I3" "I4" "I5" "I6" "I7"
## [16] "A1" "A2" "A3" "A4" "A5" "A6" "A7" "S1" "S2" "S3" "S4" "S5" "S6" "S7" "E1"
## [31] "E2" "E3" "E4" "E5" "E6" "E7" "C1" "C2" "C3" "C4" "C5" "C6" "C7"
##str(RIASEC42)
library(blcfa)
mohinh <- '
R =~ R1 + R2 + R3 + R4 + R5 + R6 + R7;
I =~ I1 + I2 + I3 + I4 + I5 + I6 + I7;
A =~ A1 + A2 + A3 + A4 + A5 + A6 + A7;
S =~ S1 + S2 + S3 + S4 + S5 + S6 + S7;
E =~ E1 + E2 + E3 + E4 + E5 + E6 + E7;
C =~ C1 + C2 + C3 + C4 + C5 + C6 + C7;
# Thêm các tham số cân nhắc cho mô hình BCFA
R ~~ I; R ~~ A; R ~~ S; R ~~ E; R ~~ C;
I ~~ A; I ~~ S; I ~~ E; I ~~ C;
A ~~ S; A ~~ E; A ~~ C;
S ~~ E; S ~~ C;
E ~~ C;
'
fit <- bcfa(model = mohinh, data = RIASEC42,
cp = "srs", n.chains = 1, burnin = 1000, sample = 5000, mcmcfile = TRUE,
mcmcextra = list(), inits = "simple",
convergence = "flexible", seed = 12345,
target = "stan", wiggle.sd = 0.1, prisamp = FALSE, jags.ic = FALSE)
##
## SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.008681 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 86.81 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 1: Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 1: Iteration: 1001 / 6000 [ 16%] (Sampling)
## Chain 1: Iteration: 1600 / 6000 [ 26%] (Sampling)
## Chain 1: Iteration: 2200 / 6000 [ 36%] (Sampling)
## Chain 1: Iteration: 2800 / 6000 [ 46%] (Sampling)
## Chain 1: Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1: Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 1: Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 1: Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 1: Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 1: Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 167.398 seconds (Warm-up)
## Chain 1: 1501.67 seconds (Sampling)
## Chain 1: 1669.07 seconds (Total)
## Chain 1:
## Computing post-estimation metrics (including lvs if requested)...
summary(fit)
## blavaan 0.5.4 ended normally after 5000 iterations
##
## Estimator BAYES
## Optimization method MCMC
## Number of model parameters 99
##
## Number of observations 1308
##
## Statistic MargLogLik PPP
## Value -73995.458 0.000
##
## Parameter Estimates:
##
##
## Latent Variables:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## R =~
## R1 1.000
## R2 1.140 0.072 1.005 1.288 1.000 normal(0,10)
## R3 1.039 0.075 0.896 1.191 1.000 normal(0,10)
## R4 0.731 0.067 0.609 0.865 1.000 normal(0,10)
## R5 0.507 0.067 0.378 0.640 1.000 normal(0,10)
## R6 0.814 0.076 0.669 0.973 1.000 normal(0,10)
## R7 0.066 0.069 -0.071 0.202 1.000 normal(0,10)
## I =~
## I1 1.000
## I2 1.170 0.048 1.080 1.267 1.000 normal(0,10)
## I3 1.074 0.046 0.989 1.167 1.000 normal(0,10)
## I4 0.849 0.053 0.750 0.954 1.000 normal(0,10)
## I5 1.066 0.048 0.974 1.167 1.000 normal(0,10)
## I6 0.900 0.051 0.802 1.002 1.000 normal(0,10)
## I7 0.737 0.050 0.643 0.842 1.000 normal(0,10)
## A =~
## A1 1.000
## A2 2.014 0.290 1.561 2.701 1.000 normal(0,10)
## A3 0.839 0.145 0.596 1.168 1.000 normal(0,10)
## A4 2.225 0.328 1.705 2.997 1.000 normal(0,10)
## A5 1.895 0.286 1.448 2.554 1.000 normal(0,10)
## A6 1.610 0.254 1.202 2.224 1.000 normal(0,10)
## A7 2.362 0.346 1.816 3.179 1.000 normal(0,10)
## S =~
## S1 1.000
## S2 1.121 0.064 1.003 1.253 1.001 normal(0,10)
## S3 1.013 0.057 0.907 1.132 1.000 normal(0,10)
## S4 0.973 0.067 0.845 1.110 1.000 normal(0,10)
## S5 1.046 0.061 0.931 1.172 1.000 normal(0,10)
## S6 1.050 0.067 0.920 1.184 1.001 normal(0,10)
## S7 1.016 0.066 0.891 1.149 1.000 normal(0,10)
## E =~
## E1 1.000
## E2 1.165 0.075 1.028 1.327 1.000 normal(0,10)
## E3 1.051 0.077 0.911 1.209 1.000 normal(0,10)
## E4 1.281 0.086 1.122 1.457 1.000 normal(0,10)
## E5 1.250 0.086 1.094 1.429 1.000 normal(0,10)
## E6 1.257 0.084 1.103 1.431 1.000 normal(0,10)
## E7 0.805 0.073 0.666 0.954 1.000 normal(0,10)
## C =~
## C1 1.000
## C2 0.847 0.045 0.760 0.937 1.000 normal(0,10)
## C3 0.878 0.045 0.791 0.971 1.000 normal(0,10)
## C4 0.936 0.053 0.834 1.043 1.000 normal(0,10)
## C5 0.870 0.050 0.777 0.970 1.000 normal(0,10)
## C6 1.082 0.051 0.986 1.185 1.000 normal(0,10)
## C7 0.789 0.050 0.693 0.891 1.000 normal(0,10)
##
## Covariances:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## R ~~
## I 0.288 0.021 0.247 0.330 1.000 lkj_corr(1)
## A 0.097 0.014 0.069 0.126 1.000 lkj_corr(1)
## S 0.174 0.015 0.146 0.206 1.000 lkj_corr(1)
## E 0.218 0.019 0.184 0.258 1.000 lkj_corr(1)
## C 0.239 0.020 0.201 0.280 1.000 lkj_corr(1)
## I ~~
## A 0.123 0.018 0.088 0.158 1.000 lkj_corr(1)
## S 0.187 0.016 0.157 0.220 1.000 lkj_corr(1)
## E 0.257 0.021 0.218 0.299 1.000 lkj_corr(1)
## C 0.303 0.022 0.261 0.348 1.000 lkj_corr(1)
## A ~~
## S 0.088 0.014 0.061 0.117 1.000 lkj_corr(1)
## E 0.113 0.017 0.081 0.147 1.000 lkj_corr(1)
## C 0.117 0.017 0.083 0.152 1.000 lkj_corr(1)
## S ~~
## E 0.197 0.017 0.165 0.232 1.000 lkj_corr(1)
## C 0.245 0.019 0.209 0.284 1.000 lkj_corr(1)
## E ~~
## C 0.287 0.023 0.244 0.334 1.000 lkj_corr(1)
##
## Variances:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## .R1 0.568 0.025 0.520 0.620 1.000 gamma(1,.5)[sd]
## .R2 0.632 0.027 0.581 0.687 1.000 gamma(1,.5)[sd]
## .R3 0.850 0.036 0.781 0.924 1.000 gamma(1,.5)[sd]
## .R4 0.854 0.035 0.787 0.927 1.000 gamma(1,.5)[sd]
## .R5 0.986 0.040 0.909 1.067 1.000 gamma(1,.5)[sd]
## .R6 1.010 0.042 0.931 1.094 1.000 gamma(1,.5)[sd]
## .R7 1.198 0.046 1.112 1.292 1.000 gamma(1,.5)[sd]
## .I1 0.516 0.023 0.473 0.565 1.000 gamma(1,.5)[sd]
## .I2 0.363 0.019 0.328 0.402 1.000 gamma(1,.5)[sd]
## .I3 0.369 0.018 0.334 0.406 1.000 gamma(1,.5)[sd]
## .I4 0.992 0.041 0.916 1.076 1.000 gamma(1,.5)[sd]
## .I5 0.472 0.022 0.432 0.516 1.000 gamma(1,.5)[sd]
## .I6 0.782 0.033 0.720 0.846 1.000 gamma(1,.5)[sd]
## .I7 0.905 0.037 0.836 0.979 1.000 gamma(1,.5)[sd]
## .A1 1.359 0.054 1.256 1.467 1.000 gamma(1,.5)[sd]
## .A2 0.751 0.034 0.684 0.820 1.000 gamma(1,.5)[sd]
## .A3 0.811 0.032 0.750 0.877 1.000 gamma(1,.5)[sd]
## .A4 0.734 0.035 0.668 0.807 1.000 gamma(1,.5)[sd]
## .A5 1.038 0.045 0.952 1.129 1.000 gamma(1,.5)[sd]
## .A6 1.063 0.045 0.978 1.155 1.000 gamma(1,.5)[sd]
## .A7 0.595 0.031 0.536 0.658 1.000 gamma(1,.5)[sd]
## .S1 0.543 0.025 0.495 0.593 1.000 gamma(1,.5)[sd]
## .S2 0.673 0.031 0.614 0.735 1.000 gamma(1,.5)[sd]
## .S3 0.496 0.023 0.452 0.543 1.000 gamma(1,.5)[sd]
## .S4 0.827 0.036 0.761 0.900 1.000 gamma(1,.5)[sd]
## .S5 0.532 0.025 0.486 0.582 1.000 gamma(1,.5)[sd]
## .S6 0.791 0.035 0.724 0.860 1.000 gamma(1,.5)[sd]
## .S7 0.738 0.032 0.678 0.803 1.000 gamma(1,.5)[sd]
## .E1 0.876 0.036 0.807 0.949 1.000 gamma(1,.5)[sd]
## .E2 0.599 0.027 0.548 0.655 1.000 gamma(1,.5)[sd]
## .E3 0.871 0.037 0.801 0.944 1.000 gamma(1,.5)[sd]
## .E4 0.724 0.033 0.661 0.792 1.000 gamma(1,.5)[sd]
## .E5 0.754 0.033 0.692 0.820 1.000 gamma(1,.5)[sd]
## .E6 0.685 0.030 0.628 0.746 1.000 gamma(1,.5)[sd]
## .E7 1.147 0.046 1.061 1.241 1.000 gamma(1,.5)[sd]
## .C1 0.562 0.027 0.510 0.615 1.000 gamma(1,.5)[sd]
## .C2 0.644 0.028 0.590 0.702 1.000 gamma(1,.5)[sd]
## .C3 0.615 0.029 0.561 0.673 1.000 gamma(1,.5)[sd]
## .C4 0.853 0.038 0.781 0.932 1.000 gamma(1,.5)[sd]
## .C5 0.767 0.034 0.703 0.836 1.000 gamma(1,.5)[sd]
## .C6 0.649 0.030 0.593 0.708 1.000 gamma(1,.5)[sd]
## .C7 0.930 0.038 0.856 1.006 1.000 gamma(1,.5)[sd]
## R 0.242 0.025 0.195 0.293 1.000 gamma(1,.5)[sd]
## I 0.441 0.034 0.376 0.509 1.000 gamma(1,.5)[sd]
## A 0.093 0.024 0.049 0.145 1.000 gamma(1,.5)[sd]
## S 0.317 0.029 0.262 0.374 1.001 gamma(1,.5)[sd]
## E 0.284 0.032 0.225 0.351 1.000 gamma(1,.5)[sd]
## C 0.501 0.039 0.429 0.581 1.000 gamma(1,.5)[sd]
Including Plots
library(semPlot)
semPaths(fit, whatLabels = "est", what = "equality",
layout="circle2", rotation = 1, style="lisrel",
bifactor="global",edge.label.cex = 1,
groups = "latents",
pastel = T, mar = c(3, 2, 2, 2))
