devtools::load_all(".")
## ℹ Loading cSEM
library(cSEM)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# Sample.R

model_Bergami="
# Measurement models
OrgPres <~ cei1 + cei2 + cei3 + cei4 + cei5 + cei6 + cei7 + cei8 
OrgIden =~ ma1 + ma2 + ma3 + ma4 + ma5 + ma6
AffJoy =~ orgcmt1 + orgcmt2 + orgcmt3 + orgcmt7
AffLove  =~ orgcmt5 + orgcmt6 + orgcmt8
Gender<~ gender

Aff<~AffJoy+AffLove

# Structural model 
OrgIden ~ OrgPres 
Aff ~ OrgPres+OrgIden+Gender
"


out <- csem(.data = BergamiBagozzi2000,.model = model_Bergami,
                       .disattenuate = T,
                       .PLS_weight_scheme_inner = 'factorial',
                       .tolerance = 1e-5,
                       .resample_method = 'bootstrap',.R = 499)

graph_attrs = c(
  # Coloring specific nodes
  "OrgPres [style=filled, fillcolor=orange]",
  "Gender [style=filled, fillcolor=orange]",
  "OrgIden [style=filled, fillcolor=lightblue]",
  "AffLove [style=filled, fillcolor=lightblue]",
  "AffJoy [style=filled, fillcolor=lightblue]",
  
  # Adjusting global graph attributes, e.g., rank separation
  "rankdir=LR",
  "ranksep=1.5",
  "nodesep=0.5",
  "splines=line",
  
  # Setting background color
  #"bgcolor=lightgray",
  
  # Vertically aligning certain nodes (If rankdir=TB, align them horizontally instead.)
  "{rank=same; OrgPres; Gender}",
  "{rank=same; OrgIden}",
  "{rank=same; AffJoy; AffLove}"
)

plot(out)
plot(out,
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "none")
plot(out,
     .plot_correlations = "none",
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "exo")
plot(out,
     .plot_correlations = "exo",
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "convcv")
plot(out,
     .plot_correlations = "convcv",
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "convcv",
     .graph_attrs = graph_attrs)
plot(out,
     .plot_correlations = "convcv",
     .graph_attrs = graph_attrs,
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "exoind")
plot(out,
     .plot_correlations = "exoind",
     .plot_labels = FALSE)
plot(out,
     .plot_correlations = "indvcv")
plot(out,
     .plot_correlations = "indvcv",
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

# 1. Non-linear models
library(moments)

## 1.1 Second-order Construct with Interaction
rm(list=ls())
model1=
'ETA3 ~ ETA1 + ETA2 + ETA1.ETA2
ETA1 =~ Y1 + Y2 + Y3
Y1=~y11+y12
Y2=~y21+y22+y23+y24
Y3=~y31+y32+y33+y34+y35+y36+y37+y38
ETA2 =~ y4 + y5 + y6
ETA3 =~ y7 + y8 + y9
'

# prameters you can choose
gam3=c();gam4=c();gam5=c()
gam3[1]=0.2 #eta1
gam3[2]=0.1 #eta2
gam3[3]=0.3 #eta1eta2

rho12=.4

# Sigma between variables on the RHS 1. equation
Sigmaexplana1 = matrix(c(1, rho12, 0,
                         rho12, 1, 0,
                         0, 0, 1+rho12^2), ncol = 3, dimnames=list(c('ETA1','ETA2','ETA1.ETA2'),c('ETA1','ETA2','ETA1.ETA2')))
Psi1=1-t(gam3)%*%Sigmaexplana1%*%t(t(gam3))

# factor loadings
Load1=c(0.6,0.7,0.8) #eta1
Load2=c(0.5,0.8,0.8) #eta2
Load3=c(0.8,0.8,0.8) #eta3

Sigma=diag(c(Psi1,1-Load1^2,1-Load2^2,1-Load3^2))

dimnames(Sigma)=list(c('Zeta1','eps1','eps2','eps3','eps4','eps5','eps6','eps7','eps8','eps9'),
                     c('Zeta1','eps1','eps2','eps3','eps4','eps5','eps6','eps7','eps8','eps9'))

# eta1 is a second-order common factorfactor
Load1_1=c(0.8,0.6)
Load1_2=c(0.8,0.7,0.6,0.8)
Load1_3=c(0.6,0.5,0.5,0.7,0.8,0.9,0.7,0.8)

Sigmafirst=diag(c(1-Load1_1^2,1-Load1_2^2,1-Load1_3^2))
dimnames(Sigmafirst)=list(c(paste0('eps1',1:2),paste0('eps2',1:4),paste0('eps3',1:8)),
                          c(paste0('eps1',1:2),paste0('eps2',1:4),paste0('eps3',1:8)))

nObs=100000

X <- MASS::mvrnorm(nObs, mu = rep(0,3),Sigma = Sigmaexplana1, empirical = F)

eta1 <- X[,1]
eta2 <- X[,2]

datexog=MASS::mvrnorm(nObs,mu=rep(0,10),Sigma=Sigma,empirical=F)
datexogfirst=MASS::mvrnorm(nObs,mu=rep(0,14),Sigma=Sigmafirst,empirical=F)

# create scores for the first-order constructs
dat1=eta1%*%t(Load1)+datexog[,c('eps1','eps2','eps3')]
colnames(dat1)=c('y1','y2','y3')
# Create observations for the indictaors connected to first-order constructs
dat11=dat1[,'y1']%*%t(Load1_1) +datexogfirst[,paste0('eps1',1:2)]
colnames(dat11)=paste0('y1',1:2)
dat12=dat1[,'y2']%*%t(Load1_2) +datexogfirst[,paste0('eps2',1:4)]
colnames(dat12)=paste0('y2',1:4)
dat13=dat1[,'y3']%*%t(Load1_3) +datexogfirst[,paste0('eps3',1:8)]
colnames(dat13)=paste0('y3',1:8)

dat2=eta2%*%t(Load2)+datexog[,c('eps4','eps5','eps6')]
dimnames(dat2)=list(c(),c('y4','y5','y6'))

# create observations for eta3
eta3 = gam3[1]*eta1 + gam3[2]*eta2 +gam3[3]*(eta1*eta2 - rho12) +datexog[,'Zeta1']
# eta3star = gam3[1]*eta1 + gam3[2]*eta2 +gam3[3]*(eta1*eta2 - rho12)

# create observations for the indicators connected to eta5
dat3=eta3%*%t(Load3)+datexog[,c('eps7','eps8','eps9')]
dimnames(dat3)=list(c(),c('y7','y8','y9'))

dat=cbind(dat11,dat12,dat13,dat2,dat3)

out=csem(.data=dat,.model=model1,.normality=F)

plot(out,
     .title = "Second-order Construct with Interaction")
plot(out,
     .title = "Second-order Construct with Interaction",
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 1.2 Second-order Construct with Multi-group
rm(list=ls())
modelSOCMMGA='
  OrgPres1 <~ cei1 + cei2 + cei3 
  OrgPres2 <~ cei4 + cei5 + cei6
  OrgPres3 =~ cei7 + cei8 
  OrgPres<~OrgPres1+OrgPres2+OrgPres3
  OrgIden =~ ma1 + ma2 + ma3 + ma4 + ma5 + ma6
  AffJoy =~ orgcmt1 + orgcmt2 + orgcmt3 + orgcmt7
  AffLove  =~ orgcmt5 + orgcmt6 + orgcmt8
  
  # Structural model 
  OrgIden ~ OrgPres 
  AffLove ~ OrgPres+OrgIden
  AffJoy  ~ OrgPres+OrgIden'

# Estimate the model for both gender groups, i.e., gender = 1 and gender = 2, using the .id argument.
# In fact, in the background the dataset is split with regard to gender and 
# the model is estimated twice.
outBergami <- csem(.data = BergamiBagozzi2000,.model = modelSOCMMGA,
                   .disattenuate = T,
                   .PLS_weight_scheme_inner = 'factorial',
                   .tolerance = 1e-6,
                   .resample_method = 'bootstrap',
                   .id='gender',
                   .seed = 1234)

# Multiple plots are generated based on the number of groups, 
# with each group having its own plot.

graph_attrs <- c(
  "rankdir=LR",
  "nodesep=0.5",
  "{rank=same; OrgPres1; OrgPres2; OrgPres3}",
  "{rank=same; AffLove;AffJoy}"
)

plot(outBergami,
     .title = "Second-order Construct with Multi-group")
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .title = "Second-order Construct with Multi-group",
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .plot_structural_model_only = TRUE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
## 1.3 Multi-group with Interaction
rm(list=ls())
model_Bergami_int="
# Measurement models
OrgPres <~ cei1 + cei2 + cei3 + cei4 + cei5 + cei6 + cei7 + cei8 
OrgIden =~ ma1 + ma2 + ma3 + ma4 + ma5 + ma6
AffJoy =~ orgcmt1 + orgcmt2 + orgcmt3 + orgcmt7
AffLove  =~ orgcmt5 + orgcmt6 + orgcmt8

# Structural model 
OrgIden ~ OrgPres 
AffLove ~ OrgPres+OrgIden+OrgPres.OrgIden
AffJoy  ~ OrgPres+OrgIden"

# Estimate the model for both gender groups, i.e., gender = 1 and gender = 2, using the .id argument.
# In fact, in the background the dataset is split with regard to gender and 
# the model is estimated twice.
outBergamiInt <- csem(.data = BergamiBagozzi2000,.model = model_Bergami_int,
                      .disattenuate = T,
                      .PLS_weight_scheme_inner = 'factorial',
                      .tolerance = 1e-6,
                      .resample_method = 'bootstrap',
                      .id='gender')

# Multiple plots are generated.

plot(outBergamiInt,
     .title = "Multi-group with Interaction")
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergamiInt,
     .title = "Multi-group with Interaction",
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergamiInt,
     .plot_structural_model_only = TRUE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergamiInt,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 1.4. Second-order Construct + Interaction + Multi-group
rm(list=ls())
modelSOCMMGA2='
  OrgPres1 <~ cei1 + cei2 + cei3 
  OrgPres2 <~ cei4 + cei5 + cei6
  OrgPres3 =~ cei7 + cei8 
  OrgPres<~OrgPres1+OrgPres2+OrgPres3
  OrgIden =~ ma1 + ma2 + ma3 + ma4 + ma5 + ma6
  AffJoy =~ orgcmt1 + orgcmt2 + orgcmt3 + orgcmt7
  AffLove  =~ orgcmt5 + orgcmt6 + orgcmt8
  
  # Structural model 
  OrgIden ~ OrgPres 
  AffLove ~ OrgPres+OrgIden
  AffJoy  ~ OrgPres+OrgIden+OrgPres.OrgIden'

outBergami2 <- csem(.data = BergamiBagozzi2000,.model = modelSOCMMGA2,
                    .disattenuate = T,
                    .PLS_weight_scheme_inner = 'factorial',
                    .tolerance = 1e-6,
                    .resample_method = 'bootstrap',
                    .id='gender',
                    .seed = 1234)

plot(outBergami2,
     .title = "Second-order Construct + Interaction + Multi-group")
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami2,
     .title = "Second-order Construct + Interaction + Multi-group",
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot_obj <- plot(outBergami2,
     .title = "Second-order Construct + Interaction + Multi-group",
     .plot_labels = FALSE)

plot(outBergami2,
     .plot_structural_model_only = TRUE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami2,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
# Testing save function

savePlot(plot_obj, .filename = "Multigroup Model with Interaction.pdf")
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_1.pdf
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_2.pdf
savePlot(plot_obj, .filename = "Multigroup Model with Interaction.png")
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_1.png
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_2.png
savePlot(plot_obj, .filename = "Multigroup Model with Interaction.dot")
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_1.dot
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_2.dot
savePlot(plot_obj, .filename = "Multigroup Model with Interaction.svg")
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_1.svg
## Plot saved to /Users/maxm1/cSEM/Multigroup Model with Interaction_2.svg
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 2.CCA.R
rm(list=ls())
model_IT_Flex="
# Composite models
ITComp  <~ ITCOMP1 + ITCOMP2 + ITCOMP3 + ITCOMP4
Modul   <~ MOD1 + MOD2 + MOD3 + MOD4
ITConn  <~ ITCONN1 + ITCONN2 + ITCONN3 + ITCONN4
ITPers  <~ ITPSF1 + ITPSF2 + ITPSF3 + ITPSF4

# Saturated structural model
ITPers ~ ITComp + Modul + ITConn
Modul  ~ ITComp + ITConn 
ITConn ~ ITComp 
"

# Estimate the model using PLS:
csem_results <- csem(.data = ITFlex, 
                     .model = model_IT_Flex,
                     .resample_method = "bootstrap",
                     .dominant_indicators = c("ITComp" = "ITCOMP1","ITConn"="ITCONN1",
                                              "Modul"="MOD1","ITPers"="ITPSF1"))

plot(csem_results,
     .title = "CCA")
plot(csem_results,
     .title = "CCA",
     .plot_labels = FALSE)
plot(csem_results,
     .plot_structural_model_only = TRUE)
plot(csem_results,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 3.Second-order constructs.R
rm(list=ls())
model_HOC='
# Measurement models
PR =~ PR1 + PR2 + PR3
IM =~ IM1 + IM2 + IM3
DI =~ DI1 + DI2 + DI3
AD =~ AD1 + AD2 + AD3
DL =~ DL1 + DL2 + DL3

# Dimensions forming the second-order construct, i.e., dimensions of brand equity
AA =~ AA1 + AA2 + AA3 + AA4 + AA5 + AA6
LO =~ LO1 + LO3
QL =~ QL1 + QL2 + QL3 + QL4 + QL5 + QL6

# SOC emergent variable (brand equity)
BE <~ QL + LO + AA

# Structural model
BE~ PR + IM + DI + AD + DL
'

# Estimate the model. 
# Internally the three-stage approach proposed by Van Riel et al. (2017) is applied.
out <- csem(.data = Yooetal2000, .model = model_HOC,
            .PLS_weight_scheme_inner = 'factorial',
            .tolerance = 1e-06,
            .resample_method = 'bootstrap'
)

plot(out,
     .title = "Second-order constructs")
plot(out,
     .title = "Second-order constructs",
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 4.Multi-group analysis.R
rm(list=ls())
model_Bergami="
# Measurement models
OrgPres <~ cei1 + cei2 + cei3 + cei4 + cei5 + cei6 + cei7 + cei8 
OrgIden =~ ma1 + ma2 + ma3 + ma4 + ma5 + ma6
AffJoy =~ orgcmt1 + orgcmt2 + orgcmt3 + orgcmt7
AffLove  =~ orgcmt5 + orgcmt6 + orgcmt8

# Structural model 
OrgIden ~ OrgPres 
AffLove ~ OrgPres+OrgIden
AffJoy  ~ OrgPres+OrgIden"

outBergami_all <- csem(.data = BergamiBagozzi2000,.model = model_Bergami,
                       .disattenuate = T,
                       .PLS_weight_scheme_inner = 'factorial',
                       .tolerance = 1e-6,
                       .resample_method = 'bootstrap')

plot(outBergami_all)
plot(outBergami_all,
     .plot_labels = FALSE)
outBergami <- csem(.data = BergamiBagozzi2000,.model = model_Bergami,
                   .disattenuate = T,
                   .PLS_weight_scheme_inner = 'factorial',
                   .tolerance = 1e-6,
                   .resample_method = 'bootstrap',
                   .id='gender')

# Automatically generate multi-plots (suffix is `Group_id`)
plot(outBergami,
     .title = "Multi-group analysis")
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .title = "Multi-group analysis",
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .plot_structural_model_only = TRUE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
plot(outBergami,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
## $`1`
## 
## $`2`
## 
## attr(,"class")
## [1] "cSEMPlot_multi" "list"
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 5.Moderation.R
rm(list=ls())
model_Int <- "
# Measurement models
# Product involvement
INV =~ INV1 + INV2 + INV3 +INV4
# Satisfaction
SAT =~ SAT1 + SAT2 + SAT3
# Switching intention
INT =~ INT1 + INT2

# Structrual model containing an interaction term.
INT ~ INV + SAT + INV.SAT
"
out <- csem(.data = Switching, .model = model_Int,
            # ADANCO settings
            .PLS_weight_scheme_inner = 'factorial',
            .tolerance = 1e-06,
            .resample_method = 'bootstrap'
)

plot(out,
     .title = "Moderation")
plot(out,
     .title = "Moderation",
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 6.Mediation.R
rm(list=ls())
model_Med='
# Measurement models
Trust =~ trust1 + trust2
PrCon =~ privcon1 + privcon2 + privcon3 + privcon4
Risk =~ risk1 + risk2 + risk3
Int =~ intent1 + intent2

# Structrual model
Int ~ Trust + PrCon + Risk
Risk ~ Trust + PrCon
Trust ~ PrCon
'

# Perform estimation
out <- csem(.data = LancelotMiltgenetal2016, .model = model_Med,
            # To reproduce the ADANCO results
            .PLS_weight_scheme_inner = 'factorial',
            .tolerance = 1e-06,
            .resample_method = 'bootstrap'
)

plot(out,
     .title = "Mediation 1")
plot(out,
     .title = "Mediation 1",
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
# To obtain inference for that indirect effect we need to write small functions.
# The function requires an .object as input.
# Those functions can be passed to csem.

IndirectEffectViaTrust <- function(.object){
  .object$Estimates$Path_estimates['Trust','PrCon']*.object$Estimates$Path_estimates['Int','Trust']
}

# Trust -> Intention
IndirectEffectViaRisk <- function(.object){
  .object$Estimates$Path_estimates['Risk','PrCon']*.object$Estimates$Path_estimates['Int','Risk']
}

# Provide the user-written function to cSEM
# Consequently, we can obtain inference about the indirect effects
outuser <- csem(.data = LancelotMiltgenetal2016, .model = model_Med,
                # To reproduce the ADANCO results
                .PLS_weight_scheme_inner = 'factorial',
                .tolerance = 1e-06,
                .resample_method = 'bootstrap',
                .user_funs = list('IE1'=IndirectEffectViaTrust,'IE2'=IndirectEffectViaRisk)
)

plot(outuser,
     .title = "Mediation 2")
plot(outuser,
     .title = "Mediation 2",
     .plot_labels = FALSE)
plot(outuser,
     .plot_structural_model_only = TRUE)
plot(outuser,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 7.OrdPLS.R
rm(list=ls())
model_Med='
# Measurement models
Trust =~ trust1 + trust2
PrCon =~ privcon1 + privcon2 + privcon3 + privcon4
Risk =~ risk1 + risk2 + risk3
Int =~ intent1 + intent2

# Structrual model
Int ~ Trust + PrCon + Risk
Risk ~ Trust + PrCon
Trust ~PrCon
'
# To let cSEM know which variables are ordinal, we need specify them as such. 
# In this case, we declare the indicators of privacy concerns as ordinal
dataord=LancelotMiltgenetal2016
for(i in c('privcon1','privcon2', 'privcon3','privcon4')){
  dataord[,i]=as.ordered(LancelotMiltgenetal2016[,i])
}

# Check the first indicator of privacy concerns
dataord$privcon1
##    [1] 5 5 5 5 5 4 4 5 5 4 4 4 1 5 4 4 5 4 3 3 3 4 5 4 4 3 5 4 2 3 5 4 5 5 5 3 3
##   [38] 5 4 4 1 5 2 4 5 3 3 5 5 4 5 5 5 5 5 3 5 4 4 4 3 5 3 4 4 5 3 5 5 5 3 4 3 2
##   [75] 3 5 1 4 5 5 5 4 4 5 4 3 5 5 4 5 4 5 5 5 4 5 5 3 5 4 3 4 4 5 5 5 5 4 4 4 5
##  [112] 5 3 5 4 5 5 5 5 4 4 4 4 4 4 4 5 4 4 4 5 4 5 5 4 3 3 3 4 5 4 5 4 5 5 5 5 4
##  [149] 4 4 4 4 2 5 5 5 5 5 2 2 4 4 3 4 5 5 4 4 5 4 4 5 4 3 5 4 5 5 4 5 5 3 4 4 5
##  [186] 4 5 5 5 4 5 4 3 5 5 5 5 3 3 5 3 4 4 3 4 2 4 5 5 4 4 4 4 5 5 5 3 5 4 4 4 4
##  [223] 5 5 3 5 5 3 3 5 5 4 4 5 5 5 3 5 4 5 5 5 2 5 5 5 5 3 5 5 4 2 5 5 2 5 5 4 5
##  [260] 5 5 5 5 4 5 4 4 5 4 1 4 5 5 4 3 4 5 4 4 5 3 5 4 5 4 1 5 3 5 5 4 4 5 4 4 3
##  [297] 4 5 4 5 5 4 4 5 3 5 5 4 5 4 3 3 5 5 4 5 4 5 5 5 3 5 3 4 5 2 3 2 4 3 4 3 5
##  [334] 5 5 2 5 5 5 4 5 5 4 3 5 5 2 4 5 5 5 4 5 4 5 3 3 5 5 5 5 4 2 4 4 3 4 4 4 5
##  [371] 3 3 3 4 4 3 5 5 4 3 5 5 3 4 5 5 5 3 3 5 2 4 4 3 3 5 5 4 5 4 4 5 5 5 2 5 5
##  [408] 4 4 4 5 4 4 4 2 2 3 5 5 3 4 2 4 2 5 3 5 5 5 5 5 4 4 5 4 5 5 4 4 5 4 1 3 3
##  [445] 5 5 5 5 5 5 4 5 5 5 5 5 3 5 5 4 3 3 3 4 5 5 5 3 1 4 5 4 3 4 5 4 5 4 5 5 3
##  [482] 2 5 3 4 2 5 3 5 5 3 5 4 4 3 4 4 5 5 5 5 5 5 4 5 5 4 5 5 4 4 4 5 5 4 2 5 5
##  [519] 3 4 5 5 4 5 5 5 4 5 5 2 5 3 3 5 5 5 4 5 5 5 5 4 5 5 5 5 5 5 5 5 4 5 5 4 3
##  [556] 5 5 4 4 5 5 5 4 5 4 4 3 4 5 5 3 5 4 5 5 3 3 4 3 3 3 3 5 4 5 4 5 4 4 4 4 4
##  [593] 5 5 5 5 4 5 5 5 4 3 5 5 5 5 5 4 1 5 5 5 4 4 5 5 5 1 5 4 3 4 5 5 5 3 5 4 5
##  [630] 2 3 4 5 5 5 4 2 5 5 5 4 5 5 5 5 5 5 5 5 3 5 4 4 3 4 5 4 3 4 5 5 5 5 3 5 4
##  [667] 5 3 5 4 4 5 5 4 5 4 5 5 2 5 5 4 5 5 5 4 5 4 5 5 5 5 5 5 5 5 5 5 3 5 3 5 5
##  [704] 5 5 4 5 2 5 4 5 5 4 4 5 5 3 5 1 5 2 5 4 5 5 4 5 5 5 3 3 5 5 5 3 5 5 5 4 2
##  [741] 1 4 5 4 5 4 5 5 5 3 5 5 5 5 5 5 5 2 3 5 3 5 5 5 4 4 5 4 4 5 3 5 3 5 4 5 5
##  [778] 5 5 1 4 5 3 5 3 4 5 5 5 5 5 5 4 5 5 5 5 5 5 5 4 5 3 5 5 5 5 5 5 5 3 5 5 5
##  [815] 5 4 5 3 5 3 5 5 5 5 5 5 5 5 3 5 4 5 5 5 2 5 5 5 5 5 5 5 3 5 4 5 2 3 1 5 5
##  [852] 4 5 4 4 4 2 4 5 5 5 4 5 5 5 5 5 5 5 5 3 3 5 5 5 4 2 3 5 5 5 5 5 5 5 5 5 3
##  [889] 5 5 5 5 3 5 3 5 2 5 5 5 5 5 5 5 4 5 5 5 5 4 3 3 2 5 3 2 5 5 4 3 4 5 5 4 3
##  [926] 4 5 3 5 4 5 4 5 5 5 5 4 5 5 5 5 5 5 2 4 4 5 4 5 4 5 1 4 4 3 4 2 4 4 4 4 3
##  [963] 4 4 2 2 5 4 2 1 1 4 1 4 5 5 5 3 3 3 5 4 4 5 4 5 3 5 3 5 5 4 2 5 4 5 4 5 4
## [1000] 4 4 3 4 4 5 3 4 5 3 1 4 4 5 5 5 5 5 5 4 5 5 5 4 3 5 5 5 5 5 3 5 5 5 5 3 2
## [1037] 5 5 4 2 5 3 5 5 5 5 4 5 4 5 4 5 4 5 4 5 5 4 4 2 5 5 5 4 4 4 3 4 4 3 3 5 5
## [1074] 5 5 4 5 5 5 4 1 5 5 5 4 5 5 4 5 4
## Levels: 1 < 2 < 3 < 4 < 5
outord <- csem(.data = dataord, .model = model_Med,
               # To reproduce the ADANCO results
               .PLS_weight_scheme_inner = 'factorial',
               .tolerance = 1e-06,
               .resample_method = 'bootstrap',
               .eval_plan = 'multicore' #Linux/MacOS: 'multicore'; Win: "multisession"
)
## Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
## processing ('multicore') is not supported when running R from RStudio because
## it is considered unstable. For more details, how to control forked processing
## or not, and how to silence this warning in future R sessions, see
## ?parallelly::supportsMulticore
plot(outord,
     .title = "OrdPLS")
plot(outord,
     .title = "OrdPLS",
     .plot_labels = FALSE)
plot(outord,
     .plot_structural_model_only = TRUE)
plot(outord,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 8.IPA.R
rm(list=ls())
# The original indicators ranged from 1-6 (6 point scale). They have been transformed to 0 - 100:
# xnew = (xold - 1) * 100/(6-1)
head(SQ)
##   a01 a02 a03 a04 a05 a06 a07 a08 a09 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19
## 1  80  80  80  80  80  80  80  80 100  80  80 100 100 100 100 100 100  80  80
## 2  80  80  80  80 100  80  80  80 100 100 100 100  80  80  80  80  60  80  80
## 3  80  80 100  80 100  80  80 100 100  80  80  80  80  80  80 100  80  80 100
## 4  80  80  80 100 100 100  80 100 100 100 100 100 100 100 100  80  80 100  80
## 5  60  60  40  60 100  80 100 100 100 100 100 100  80 100 100 100 100 100 100
## 6 100  80  80  80  80  80  80  80  80  80  80  80  80  80  80  80  80  80  80
##   a20 a21 a22 sat
## 1  80  80  80  80
## 2  80  60  60  80
## 3 100  80  80  80
## 4  80  60  80  80
## 5 100 100 100  80
## 6  80  60  80  80
# Specify the model
model<-'
# Specify composite models
Tangibles <~ a01+a02+a03+a04
Reliability<~ a05+a06+a07+a08+a09
Responsiveness <~ a10+a11+a12+a13
Assurance<~a14+a15+a16+a17
Empathy<~a18+a19+a20+a21+a22
Sat <~ sat

# Structural model
Sat~Tangibles+Reliability+Responsiveness+Assurance+Empathy
'

# Estimate model
out <- csem(.data = SQ,.model = model,
            .PLS_weight_scheme_inner = 'factorial',
            .tolerance = 1e-06)

plot(out,
     .title = "IPA")
plot(out,
     .title = "IPA",
     .plot_labels = FALSE)
plot(out,
     .plot_structural_model_only = TRUE)
plot(out,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)
# Perform Importance-Performance Analysis
outIPA <- doIPMA(out)

# Illustrate IPMA

# On construct level
plot(x = outIPA,.dependent = 'Sat',.level = 'construct',
     .attributes = c("Tangibles","Reliability","Responsiveness", "Assurance","Empathy"))

# On indicator level
plot(x = outIPA , .dependent = 'Sat', .level = 'indicator',
     .attributes = c('a01', 'a02', 'a03', 'a04', 'a05', 'a06',
                     'a07', 'a08', 'a09', 'a10', 'a11', 'a12',
                     'a13', 'a14', 'a15', 'a16', 'a17', 'a18',
                     'a19', 'a20', 'a21', 'a22'))

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# 9.Sigma_Summers_composites
rm(list=ls())
model <- "
ETA1 ~ ETA2 + XI1 + XI2
ETA2 ~ ETA1 + XI3 +XI4
ETA1 ~~ ETA2
XI1 <~ x1 + x2 + x3
XI2 <~ x4 + x5 + x6
XI3 <~ x7 + x8 + x9
XI4 <~ x10 + x11 + x12
ETA1 <~ y1 + y2 + y3
ETA2 <~ y4 + y5 + y6
"
## Generate data
summers_dat <- MASS::mvrnorm(n = 300, mu = rep(0, 18), Sigma = Sigma_Summers_composites, empirical = TRUE)

## Estimate
res <- csem(.data = summers_dat, .model = model) # inconsistent
plot(res,
     .title = "Sigma_Summers_composites 1")
plot(res,
     .title = "Sigma_Summers_composites 1",
     .plot_labels = FALSE)
# 2SLS
res_2SLS <- csem(.data = summers_dat, .model = model, .approach_paths = "2SLS",
                 .instruments = list(ETA1 = c('XI1','XI2','XI3','XI4'),
                                     ETA2 = c('XI1','XI2','XI3','XI4')))

plot(res_2SLS,
     .title = "Sigma_Summers_composites 2")
plot(res_2SLS,
     .title = "Sigma_Summers_composites 2",
     .plot_labels = FALSE)
plot(res_2SLS,
     .plot_structural_model_only = TRUE)
plot(res_2SLS,
     .plot_structural_model_only = TRUE,
     .plot_labels = FALSE)