#This R package performs PLS-SEM mode A and the transfromed mode A-->B
#The package contains three functions 
#1 reg_swc()--estimate the regression coefficients represented by the solid arrows in Figure 6;
#2 weight_Ba()--transform the weights under mode A to mode B
#3 PLS_AB()--compute the weights under PLS-SEM mode A, and call the weight_Ba() function for w_A to w_Ba
#The package uses Weston & Gore Jr (2006) as an example
#It does not need raw data
# Users need to be familiar with R in order to modify the package for other
#models

#----------------------------------------------------
#Estimate the regression coefficients
#of the structural model represented by the solid arrow in Figure 6
#using the sample covariance matrix swc.cov of the weighted composites 
reg_swc=function(N,swc.cov){
	n=N-1
	hg11=swc.cov[2,1]/swc.cov[1,1] #\hat\gamma_11
	hsigma_res2=n*(swc.cov[2,2]-swc.cov[2,1]*swc.cov[1,2]/swc.cov[1,1])/(N-2) #residual variance
	Var_g11=hsigma_res2/(n*swc.cov[1,1]) #variance of hat\gamma_11
	SE_g11=sqrt(Var_g11)
	z_g11=hg11/SE_g11
	gamma_11=cbind(hg11,SE_g11,z_g11)

	hgb21=solve(swc.cov[1:2,1:2])%*%swc.cov[1:2,3] #(\hat\gamma_21,\hat\beta_21)
	hsigma_res3=n*(swc.cov[3,3]-swc.cov[3,1:2]%*%solve(swc.cov[1:2,1:2])%*%swc.cov[1:2,3])/(N-3)
	Var_gb21=c(hsigma_res3)*solve(n*swc.cov[1:2,1:2])
	SE_gb21=sqrt(diag(Var_gb21))
	z_gb21=hgb21/SE_gb21
	gammabeta_21=cbind(hgb21,SE_gb21,z_gb21)

	hb32=swc.cov[4,3]/swc.cov[3,3] #\hat\beta_32
	hsigma_res4=n*(swc.cov[4,4]-swc.cov[4,3]*swc.cov[3,4]/swc.cov[3,3])/(N-2)
	Var_b32=hsigma_res4/(n*swc.cov[3,3])
	SE_b32=sqrt(Var_b32)
	z_b32=hb32/SE_b32
	beta_32=cbind(hb32,SE_b32,z_b32)
	estimate=rbind(gamma_11,gammabeta_21,beta_32)
	rownames(estimate)=c("gamma_11", "gamma_21", "beta_21", "beta_32")
	colnames(estimate)=c("est.", "SE", "z")
	return(estimate)
}

#----------------------------------------------------
#Transform weight of model A to B according to a one-factor 
#model Sigma=h*w_A*w_A'+Psi for a single block;
#w_A is the input vector of weights obtained under mode A 
#for a block of indicators
#S is the input sample covariance matrix for the block
#corresponding to w_A
#w_Ba is the transformed weights from mode A to mode B
weight_Ba=function(w_A,S) {
p=nrow(S);
w_At=t(w_A);
c2=w_At%*%S%*%w_A; 
#c2 should be equal to 1.0 if the corresponding 
#composite is standardized, i.e. Var(w_A'x)=1;
c2=c2[1,1];
d_s=diag(S);
t1=c2-sum(w_A*w_A*d_s); #elementwise multiplication;
t2=sum(w_A^2)*sum(w_A^2)-sum(w_A^4);
h=t1/t2; 
load2=h*(w_A*w_A); #vector of squared loadings; 
psi=d_s-load2; #error variance;
for (j in 1:p) {
   psi_j=psi[j];
   if (psi_j<0) {
      psi[j]=.05; #negative error variance is replaced by .05;
      cat("Heywood case: variable", j, "\n") #output the warning;
      cat("has error variance", psi_j, "\n")
      cat("which was changed to .05 by the program", "\n")
   }
}
w_s=w_A/psi; #The transformed weight is proportional to w_a/psi;
w_st=t(w_s);
wSw=w_st%*%S%*%w_s;
w_Ba=w_s/c(sqrt(wSw)); #The weights in w_Ba satisfy Var(w_Ba'x)=1;
return(w_Ba)
}; #end of the function;

#----------------------------------------------------
#Perform PLS-SEM with a sample covariance covariance matrix, 
#which is assuemd (divided by N-1) to be unbiased estimate of 
#the poipulation matrix Sigma
#The sample covariance matrix (re-orded) scov is from Weston & Gore Jr (2006)
#The structural model is represented by the solid arrow in Figure 6

PLS_ABa=function(N,scov){
n=N-1
p=ncol(scov)
scov=n*scov/N; #rescale the sample covariance matrix for the convenience of LS regression

#put the sample covariance matrix in blocks
p_x=3;
p_y1=3;
p_y2=3;
p_y3=3;
s_xx=scov[1:3,1:3];s_xy1=scov[1:3,4:6];s_xy2=scov[1:3,7:9];s_xy3=scov[1:3,10:12];
s_y1x=t(s_xy1);      s_y11=scov[4:6,4:6];s_y12=scov[4:6,7:9];s_y13=scov[4:6,10:12];
s_y2x=t(s_xy2);      s_y21=t(s_y12);       s_y22=scov[7:9,7:9];s_y23=scov[7:9,10:12];
s_y3x=t(s_xy3);      s_y31=t(s_y13);       s_y32=t(s_y23);       s_y33=scov[10:12,10:12];

#starting value of the weights;
w_x=matrix(1,p_x,1)
c2_x=t(w_x)%*%s_xx%*%w_x;
w_x=w_x/c(sqrt(c2_x));

w_y1=matrix(1,p_y1,1)
c2_y1=t(w_y1)%*%s_y11%*%w_y1;
w_y1=w_y1/c(sqrt(c2_y1));

w_y2=matrix(1,p_y2,1)
c2_y2=t(w_y2)%*%s_y22%*%w_y2;
w_y2=w_y2/c(sqrt(c2_y2));

w_y3=matrix(1,p_y3,1)
c2_y3=t(w_y3)%*%s_y33%*%w_y3;
w_y3=w_y3/c(sqrt(c2_y3));

#the sample correlation for environment variables;
r_xy1=t(w_x)%*%s_xy1%*%w_y1;
r_xy2=t(w_x)%*%s_xy2%*%w_y2;
r_y12=t(w_y1)%*%s_y12%*%w_y2;
r_y23=t(w_y2)%*%s_y23%*%w_y3;

w_0=rbind(w_x,w_y1,w_y2,w_y3) #put the weigthts in vectors for update
#Iteratively computing mode A;
dw=1; ep=.00001; #The criterion for convergence;
for (j in 1:50){
#do j=1 to 50 while (max(abs(dw))>ep);
	#ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2; #environment variable of \xi in notation
	w_x=c(sign(r_xy1))*s_xy1%*%w_y1+c(sign(r_xy2))*s_xy2%*%w_y2;
	#w_x is the vector of regression coefficients of ev_x-->x;
	c2_x=t(w_x)%*%s_xx%*%w_x;
	w_x=c(sign(w_x[1]))*w_x/c(sqrt(c2_x)); 
	#standardize the weights so that the variance of the composite(xi_x) =1;
	
	#ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2; #environment variable of \eta_y1 in notation
	w_y1=c(sign(r_xy1))*(s_y1x)%*%w_x+c(sign(r_y12))*s_y12%*%w_y2; 
	#w_y1 is the vector of regression coefficients of ev_y1-->y1;
	c2_y1=t(w_y1)%*%s_y11%*%w_y1;
	w_y1=c(sign(w_y1[1]))*w_y1/c(sqrt(c2_y1)); #standardize the weights so that var(eta_y1)=1;
	
	#ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3; #environment variable of \eta_y2 in notation
	w_y2=c(sign(r_xy2))*s_y2x%*%w_x+c(sign(r_y12))*s_y21%*%w_y1+c(sign(r_y23))*s_y23%*%w_y3; 
	#w_y2 is the vector of regression coefficients of ev_y2-->y2;	
	c2_y2=t(w_y2)%*%s_y22%*%w_y2;
	w_y2=c(sign(w_y2[1]))*w_y2/c(sqrt(c2_y2)); 

	#ev_y3=eta_y2; #environment variable of \eta_y3 in notation
	w_y3=c(sign(r_y23))*s_y32%*%w_y2;
	#w_y3 is the vector of regression coefficients of ev_y3-->y3;
	c2_y3=t(w_y3)%*%s_y33%*%w_y3;
	w_y3=c(sign(w_y3[1]))*w_y3/c(sqrt(c2_y3));


	#update the correlation for environment variables;
	r_xy1=t(w_x)%*%s_xy1%*%w_y1;
	r_xy2=t(w_x)%*%s_xy2%*%w_y2;
	r_y12=t(w_y1)%*%s_y12%*%w_y2;
	r_y23=t(w_y2)%*%s_y23%*%w_y3;

	w_j=rbind(w_x,w_y1,w_y2,w_y3);
	dw=max(abs(w_j-w_0));
	w_0=w_j; #update the weights
	if (dw<ep){break}
} #end of loop j
if (j<=49){diverg=0}
if (j>=50){diverg=1}

w_A=matrix(0,p,4) #put the estimated weights in a matrix
w_A[1:p_x,1]=w_x
w_A[(p_x+1):(p_x+p_y1),2]=w_y1
w_A[(p_x+p_y1+1):(p_x+p_y1+p_y2),3]=w_y2
w_A[(p_x+p_y1+p_y2+1):p,4]=w_y3

#tranform weights w_A into w_Ba and put the transformed weights in a matrix
w_Ba=matrix(0,p,4)
w_bx=weight_Ba(w_x,s_xx);
w_by1=weight_Ba(w_y1,s_y11);
w_by2=weight_Ba(w_y2,s_y22);
w_by3=weight_Ba(w_y3,s_y33);
w_Ba[1:p_x,1]=w_bx
w_Ba[(p_x+1):(p_x+p_y1),2]=w_by1
w_Ba[(p_x+p_y1+1):(p_x+p_y1+p_y2),3]=w_by2
w_Ba[(p_x+p_y1+p_y2+1):p,4]=w_by3

estimate_AB=list(w_A,w_Ba,diverg)
#return(estimate_AB)

if (diverg==1){
	print("PLS-SEM mode A does not converge with 50 iterations")
}else{	
	w_A=estimate_AB[[1]]
	scA.cov=t(w_A)%*%scov%*%w_A
	print("PLS-SEM mode A")
	reg_plsA=reg_swc(N,scA.cov)
	print(reg_plsA)

	w_Ba=estimate_AB[[2]]
	scBa.cov=t(w_Ba)%*%scov%*%w_Ba
	print("PLS-SEM mode Ba")
	reg_plsBa=reg_swc(N,scBa.cov)
	print(reg_plsBa)
}


} #end of the function;