#This is the R code for path analysis using Bartlett-factor-scores 
#and equally weighted composites, with the sample covariance matrix 
#in Table 2 of Weston and Gore Jr (2006) as an example
#The order of the variables is rearranged from that in Weston and Gore Jr

source("Table7source-code.R")

#scov is the sample covariance matrix from Weston & Gore Jr
scov=scan()
7.511  5.505  5.773  3.413  2.926  2.852  0.243  0.299  0.177  3.260  3.592  2.279
5.505  7.220  5.716  2.518  2.840  2.910  0.281  0.299  0.296  3.070  3.665  2.394
5.773  5.716  6.958  2.630  2.671  2.618  0.256  0.270  0.209  2.596  2.817  2.553
3.413  2.518  2.630  4.590  3.626  3.668  0.260  0.288  0.230  2.976  2.931  2.402
2.926  2.840  2.671  3.626  3.860  3.597  0.225  0.278  0.239  2.743  2.867  2.418
2.852  2.910  2.618  3.668  3.597  4.237  0.260  0.286  0.280  2.764  2.821  2.410
0.243  0.281  0.256  0.260  0.225  0.260  0.072  0.037  0.039  0.257  0.224  0.236
0.299  0.299  0.270  0.288  0.278  0.286  0.037  0.064  0.038  0.336  0.297  0.269
0.177  0.296  0.209  0.230  0.239  0.280  0.039  0.038  0.078  0.280  0.239  0.198
3.260  3.070  2.596  2.976  2.743  2.764  0.257  0.336  0.280  4.811  3.909  3.179
3.592  3.665  2.817  2.931  2.867  2.821  0.224  0.297  0.239  3.909  5.074  3.041
2.279  2.394  2.553  2.402  2.418  2.410  0.236  0.269  0.198  3.179  3.041  3.503

scov=matrix(data = scov, nrow = 12, ncol = 12, byrow = TRUE) #make scov a matrix
N=403

#hlamb contains the estimated factor loadings for the model represented by the
#solid arrows in Figure 6 of Deng & Yuan (2022). The 1.0000 are fixed loadings 
#for scaling the endogenous constructs. The estimates can be obtained from
#EQS or another software where scov is fitted by the model using the 
#normal-distribution-based maximum likelihood.
hlamb=c(
 2.38134,
 2.36480,
 2.40249,
 1.00000,
 0.97641,
 0.99313,
 1.00000,
 1.14352,
 1.01117,
 1.00000,
 0.96297,
 0.79538
)

#hpsi contains the estimated error variances for the model represented 
#by the solid arrows in Figure 6 of Deng & Yuan (2022).
hpsi=c(
 1.84020,
 1.62771,
 1.18602,
 0.88226,
 0.32511,
 0.58005,
 0.04360,
 0.02686,
 0.04896,
 0.79454,
 1.34951,
 0.96207
)
lamb=matrix(0,12,4) #put the factor loadings in a 12 by 4 matrix;
lamb[1:3,1]=hlamb[1:3]
lamb[4:6,2]=hlamb[4:6]
lamb[7:9,3]=hlamb[7:9]
lamb[10:12,4]=hlamb[10:12]
lamb_t=t(lamb) #transpose of lamb

psi=matrix(0,12,12) #put the error variances in a 12 by 12 diagonal matrix
psi=diag(hpsi[1:12])

#compute the sample covariance matrix of the Bartlett-facor-scores
psi_in=solve(psi)
coef_t=solve(lamb_t%*%psi_in%*%lamb)%*%lamb_t%*%psi_in 
F.cov=coef_t%*%scov%*%t(coef_t) #the sample cov of Bartlett-factor-scores

print("Bartlett-factor-score regression")
#estimate the regression coefficients using reg_swc()
reg_bfs=reg_swc(N,F.cov) 
print(reg_bfs)

##############################################
#Regression with equally weighted composites
#the sum scores for each composites are divided by
#the sum of the loadings (which does not affect the z-statistic in the table)
c_ewc=matrix(0,12,4)
one_3=rep(1,3)
c_ewc[1:3,1]=one_3/sum(hlamb[1:3])
c_ewc[4:6,2]=one_3/sum(hlamb[4:6])
c_ewc[7:9,3]=one_3/sum(hlamb[7:9])
c_ewc[10:12,4]=one_3/sum(hlamb[10:12])
sc.cov=t(c_ewc)%*%scov%*%c_ewc

print("Equally-weighted-composite regression")
#estimate the regression coefficients using reg_swc()
reg_ewc=reg_swc(N,sc.cov)
print(reg_ewc)