*Data from LF Deng, and with robust transformation;
data depression;
infile 'D:\papers\2021\PLS-SEM\CB_SEM&PLS_SEM\webfile\healthstress_rob.dat';
input v1-v24;
run;
***------------------------------------------------------------***;
options linesize=100;
proc iml;
reset noname;
*----------------------------------------------------*;
**** The following subroutine is to perform vech(.) function;
*----------------------------------------------------*;
start vech(A,Va);
l=0;
p=nrow(A);
pstar=p*(p+1)/2;
Va=j(pstar,1,0);
do j=1 to p;
	do i=j to p;
    	l=l+1;
		Va[l]=A[i,j];
	end;
end;
finish;
*-----------------------------------------------------------------------*;
* creates duplication matrix as defined in Magnus&Neudecker*;
*-----------------------------------------------------------------------*;
start DP(p, dup);
  Dup=j(p*p,p*(p+1)/2,0); 
  count=0;
  do j=1 to p;
    do i=j to p;
      count=count+1;
      	if i=j then do;	
      		Dup[(j-1)*p+j, count]=1;
	end;
	else do;      
      		Dup[(j-1)*p+i, count]=1;
      		Dup[(i-1)*p+j, count]=1;
      	end;
    end;
  end;
finish;
***--------------------------------------------------------------------***;
**** Subroutine to evaluate sigma0 and (dsigma0/dtheta);
*Model as represented by Figure 8 of Deng & Yuan;
***--------------------------------------------------------------------***;
start sigdsig(theta0,vsig,Sigma,vdsig);
p=24;
p_x=15;
p_y=9;
ps=p*(p+1)/2;
lamb_x=j(15,3,0);
lamb_x[1:5,1]=theta0[1:5];
lamb_x[6:9,2]=theta0[6:9];
lamb_x[10:15,3]=theta0[10:15];
lamb_y=j(9,1,0);
lamb_y[1]=1.0;
lamb_y[2:9]=theta0[16:23];
gamma=theta0[24:26];
gamma_t=gamma`;
phi=i(3);
phi[2:3,1]=theta0[27:28]; phi[1,2:3]=(phi[2:3,1])`;
phi[3,2]=theta0[29]; phi[2,3]=phi[3,2];
sig_zeta=theta0[30];

psi_x=diag(theta0[31:45]);
psi_y=diag(theta0[46:54]);
*y=lamb_y*eta+e_y;
*x=lamb_x*xi+e_x;
*eta=gamma_t*xi+zeta;
*computing Sigma based on structured parameters;
lamb_yt=lamb_y`;
lamb_xt=lamb_x`;
sigyy=lamb_y*(gamma_t*phi*gamma+sig_zeta)*lamb_yt+psi_y;
sigxy=lamb_x*phi*gamma*lamb_yt;
sigxx=lamb_x*phi*lamb_xt+psi_x;
sigyx=sigxy`;
sigma=j(p,p,0);
sigma[1:p_x,1:p_x]=sigxx;
sigma[1:p_x,(p_x+1):p]=sigxy;
sigma[(p_x+1):p,1:p_x]=sigyx;
sigma[(p_x+1):p,(p_x+1):p]=sigyy;

run vech(sigma,vsig);

***** The following is to evaluate (dsigma0/dtheta);

* dlamb_x is the derivative with lamb_x;

dlamb_x=j(ps,p_x,0);
phgLy=phi*gamma*lamb_yt;
phLx=phi*lamb_xt;
Lxph=lamb_x*phi;

tyy=j(p_y,p_y,0);
do i=1 to 5;
	tt=j(p_x,3,0);
	tt[i,1]=1;
	txx=tt*phLx+Lxph*(tt`);
	txy=tt*phgLy;
	tyx=txy`;
	ttt=(txx||txy)//(tyx||tyy);
	run vech(ttt,tttt);
	dlamb_x[,i]=tttt;
end;
do i=6 to 9;
	tt=j(p_x,3,0);
	tt[i,2]=1;
	txx=tt*phLx+Lxph*(tt`);
	txy=tt*phgLy;
	tyx=txy`;
	ttt=(txx||txy)//(tyx||tyy);
	run vech(ttt,tttt);
	dlamb_x[,i]=tttt;
end;
do i=10 to 15;
	tt=j(p_x,3,0);
	tt[i,3]=1;
	txx=tt*phLx+Lxph*(tt`);
	txy=tt*phgLy;
	tyx=txy`;
	ttt=(txx||txy)//(tyx||tyy);
	run vech(ttt,tttt);
	dlamb_x[,i]=tttt;
end;
*dlamb_y is the derivative with lamb_y;
dlamb_y=j(ps,8,0);
gphg=gamma_t*phi*gamma+sig_zeta;
gphgLy=gphg*lamb_yt;
Lygphg=gphgLy`;
gphLx=gamma_t*phi*lamb_xt;
txx=j(p_x,p_x,0);
do i=1 to 8;
	tt=j(p_y,1,0);
	tt[i+1]=1;
	tyy=tt*gphgLy+Lygphg*(tt`);
	tyx=tt*gphLx;
	txy=tyx`;
	ttt=(txx||txy)//(tyx||tyy);
	run vech(ttt,tttt);
	dlamb_y[,i]=tttt;
end;

*dgamma is the derivative with gamma;
dgamma=j(ps,3,0);
phg=phi*gamma;
gph=gamma_t*phi;
txx=j(p_x,p_x,0);
do i=1 to 3;
	tt=j(3,1,0);
	tt[i]=1;
	tyy=lamb_y*((tt`)*phg+gph*tt)*lamb_yt;
	txy=lamb_x*phi*tt*lamb_yt;
	tyx=txy`;
	ttt=(txx||txy)//(tyx||tyy);
	run vech(ttt,tttt);
	dgamma[,i]=tttt;
end;

*dphi is the derivative with phi;
dphi=j(ps,3,0);
Lyg=lamb_y*gamma_t;
gLy=gamma*lamb_yt;

tt=j(3,3,0);
tt[2,1]=1; tt[1,2]=1;
txx=lamb_x*tt*lamb_xt;
txy=lamb_x*tt*gLy;
tyx=txy`;
tyy=Lyg*tt*gLy;
ttt=(txx||txy)//(tyx||tyy);
run vech(ttt,tttt);
dphi[,1]=tttt;

tt=j(3,3,0);
tt[3,1]=1; tt[1,3]=1;
txx=lamb_x*tt*lamb_xt;
txy=lamb_x*tt*gLy;
tyx=txy`;
tyy=Lyg*tt*gLy;
ttt=(txx||txy)//(tyx||tyy);
run vech(ttt,tttt);
dphi[,2]=tttt;

tt=j(3,3,0);
tt[3,2]=1;tt[2,3]=1;
txx=lamb_x*tt*lamb_xt;
txy=lamb_x*tt*gLy;
tyx=txy`;
tyy=Lyg*tt*gLy;
ttt=(txx||txy)//(tyx||tyy);
run vech(ttt,tttt);
dphi[,3]=tttt;

*dsig_zeta is the derivative with sig_zeta;
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
tyx=txy`;
tyy=lamb_y*lamb_yt;
ttt=(txx||txy)//(tyx||tyy);
run vech(ttt,tttt);
dsig_zeta=tttt;

*dpsi_x is the derivative with psi_x;
dpsi_x=j(ps,p_x,0);
tyy=j(p_y,p_y,0);
txy=j(p_x,p_y,0);
tyx=txy`;
do j=1 to p_x;
	tt=j(p_x,p_x,0);
    tt[j,j]=1;
    txx=tt;
    ttt=(txx||txy)//(tyx||tyy);
    run vech(ttt,tttt);
    dpsi_x[,j]=tttt;
end;

*dpsi_y is the derivative with psi_y;
dpsi_y=j(ps,p_y,0);
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
tyx=txy`;

do j=1 to p_y;
	tt=j(p_y,p_y,0);
    tt[j,j]=1;
    tyy=tt;
    ttt=(txx||txy)//(tyx||tyy);
    run vech(ttt,tttt);
    dpsi_y[,j]=tttt;
end;

vdsig=dlamb_x||dlamb_y||dgamma||dphi||dsig_zeta||dpsi_x||dpsi_y;

finish;
***--------------------------------------------------------------------***;
**Subroutine for parameter estimation by minimizing the F_ml function using Fisher-scoring method**;
**coefficients of factor-scores are subsequently obtained;
*Model as represented by Figure 8 of Deng & Yuan;
***--------------------------------------------------------------------***;
start ML_sem(n,p,dup,theta0,scov,T_ml,htheta_z,coef_xt,coef_yt,n_div);
ep=.00001;
dtheta=1;
n_div=1;
do jj=1 to 300 while (max(abs(dtheta))>ep);
	run sigdsig(theta0,vsig0,Sigma0, vdsig);
  	sig_in=inv(Sigma0);
	** weight given by normal theory ;
	weight=0.5*dup`*(sig_in@sig_in)*dup;
	vdsig_t=vdsig`;
	dsw=vdsig_t*weight;
	dswds=dsw*vdsig;
	dwd_in=inv(dswds);
	cov_res=scov-sigma0;
	run vech(cov_res,vcov_res);
	dtheta=dwd_in*dsw*vcov_res;
	theta0=theta0+dtheta;
end;

***** Test start here;
if jj<300 then do;
	n_div=0;
	Ssig_in=Scov*sig_in;
	f_ml=trace(ssig_in)-log(det(ssig_in))-p;
	T_ml=(n-1)*f_ml;
	SE_hat=sqrt(vecdiag(dwd_in)/(n-1));
	*SE_hat=sqrt(vecdiag(dwd_in)/n);
	z_score=theta0/SE_hat;
	htheta_z=theta0||SE_hat||z_score;
*computation of factor-score;
	lamb_x=j(15,3,0);
	lamb_x[1:5,1]=theta0[1:5];
	lamb_x[6:9,2]=theta0[6:9];
	lamb_x[10:15,3]=theta0[10:15];
	lamb_y=j(9,1,0);
	lamb_y[1]=1.0;
	lamb_y[2:9]=theta0[16:23];
	lamb_xt=lamb_x`;
	lamb_yt=lamb_y`;
	psi_x=diag(theta0[31:45]);
	psi_y=diag(theta0[46:54]);
	psi_xin=inv(psi_x);
	psi_yin=inv(psi_y);
	coef_xt=inv(lamb_xt*psi_xin*lamb_x)*lamb_xt*psi_xin;
	coef_yt=inv(lamb_yt*psi_yin*lamb_y)*lamb_yt*psi_yin;	
end;
finish;
*--------------------------------------------;
*LS regression using composite-scores;
*--------------------------------------------;
start regress(x,y,reg_est, Rsq);
*x has three columns, y has one column;
n=nrow(x);
*print "reg n=" n;
x_t=x`;
s_xx=x_t*x; 
s_xy=x_t*y;
s_xxin=inv(s_xx);
hbeta=s_xxin*s_xy;
hsig2=ssq(y-x*hbeta)/(n-4);
cov_beta=hsig2*s_xxin; 
SE=sqrt(vecdiag(cov_beta));
z_beta=hbeta/SE;

reg_est=hbeta||SE||z_beta;

Rsq=1-ssq(y-x*hbeta)/ssq(y);
finish;

*--------------------------------------------;
*Transforming weight of model A to B using structured 
Sigma=h*w*w'+Psi;
*--------------------------------------------;
start ssw(w,S,w_ss);
p=nrow(S);
w_t=w`;
c2=w_t*S*w;
d_s=vecdiag(S);
t1=c2-sum(w#w#d_s);
t2=ssq(w)*ssq(w)-sum(w#w#w#w);
h=t1/t2; 
psi=d_s-h*w#w;
*print "p, psi=" p psi;
do j=1 to p;
	psi_j=psi[j];
	if psi_j<0 then do;
		*print "Heywood case, j=" j psi_j;
		psi[j]=.05;
	end;
end;
w_s=w/psi;
w_ss=w_s/sqrt(w_s`*S*w_s);
finish;
******************************************************************;
*Subroutine conducting PLS-SEM mode A, and B_A;
*PLS-SEM mode A based on the raw data (centralized);
*the model is as represented by Figure 8 of Deng & Yuan;
******************************************************************;
start PLS_A(x1,x2,x3,y, reg_plsA, reg_plsBA, Rsq_A, Rsq_BA);
n=nrow(x1); 
p1=ncol(x1);
p2=ncol(x2);
p3=ncol(x3);
py=ncol(y);

s_11=x1`*x1/n;
s_22=x2`*x2/n;
s_33=x3`*x3/n;
s_yy=y`*y/n;
*--------------------;
*starting value of weights;
w_1=j(p1,1,1);
c2_1=w_1`*s_11*w_1;
w_1=w_1/sqrt(c2_1);

w_2=j(p2,1,1);
c2_2=w_2`*s_22*w_2;
w_2=w_2/sqrt(c2_2);

w_3=j(p3,1,1);
c2_3=w_3`*s_33*w_3;
w_3=w_3/sqrt(c2_3);

w_y=j(py,1,1);
c2_y=w_y`*s_yy*w_y;
w_y=w_y/sqrt(c2_y);

eta_x1=x1*w_1;
eta_x2=x2*w_2;
eta_x3=x3*w_3;
eta_y=y*w_y;
w_0=w_1//w_2//w_3//w_y;
*print "-------------mode A-------------";
dw=1; ep=.00001;
do j=1 to 50 while (max(abs(dw))>ep);
	w_1=x1`*eta_y/n;
	*w_1 is the vector of regression coefficients eta_y-->x1, var(eta_y)=1;
	c2_1=w_1`*s_11*w_1;
	w_1=sign(w_1[1])*w_1/sqrt(c2_1);
	eta_x1=x1*w_1; *eta_x1 is the vector of predicted score;

	w_2=x2`*eta_y/n;
	*w_2 is the vector of regression coefficients eta_y-->x2, var(eta_y)=1;
	c2_2=w_2`*s_22*w_2;
	w_2=sign(w_2[1])*w_2/sqrt(c2_2);
	eta_x2=x2*w_2; *eta_x2 is the vector of predicted score;

	w_3=x3`*eta_y/n;
	*w_3 is the vector of regression coefficients eta_y--> x3, var(eta_y)=1;
	c2_3=w_3`*s_33*w_3;
	w_3=sign(w_3[1])*w_3/sqrt(c2_3);
	eta_x3=x3*w_3; *eta_x3 is the vector of predicted score;

	r_y1=(eta_y`*eta_x1)/n;
	r_y2=(eta_y`*eta_x2)/n;
	r_y3=(eta_y`*eta_x3)/n;

	*environmental variable of y;
	en_y=sign(r_y1)*eta_x1+sign(r_y2)*eta_x2+sign(r_y3)*eta_x3;
	
	w_y=y`*en_y/n;*w_y is the vector of regression coefficients;
	c2_y=w_y`*s_yy*w_y;
	w_y=sign(w_y[1])*w_y/sqrt(c2_y);
	eta_y=y*w_y; *eta_y is the vector of predicted score;
	*print jj w_x1 w_x2;
	w_j=w_1//w_2//w_3//w_y;
	dw=w_j-w_0;
	w_0=w_j;
end;
print "mode A iteration time=" j;
*print "w_x1,w_x2,w_x3,w_y=" w_1 w_2 w_3 w_y;
*print "w_y,c2_y=" w_y[format=f6.3] c2_y;
w_a=w_1//w_2//w_3//w_y;
*regression of eta_y on eta_x1, eta_x2, eta_x3;
x_eta=eta_x1||eta_x2||eta_x3;
y_eta=eta_y;
run regress(x_eta,y_eta,reg_plsA, Rsq_A);

*print "--------mode A structured to mode B--------";
*print "--------block-1---------";
run ssw(w_1,S_11,w_s1); 
*print "--------block-2---------";
run ssw(w_2,S_22,w_s2);
*print "--------block-3---------";
run ssw(w_3,S_33,w_s3);
*print "--------block-y---------";
run ssw(w_y,S_yy,w_sy);


w_sb=w_s1//w_s2//w_s3//w_sy;

sx_eta=(x1*w_s1)||(x2*w_s2)||(x3*w_s3);
sy_eta=y*w_sy;
run regress(sx_eta,sy_eta,reg_plsBA, Rsq_BA);
*print "---------------------------------weights---------------------------------";
*print "mode A(w_x1,w_x2,w_x3,w_y)--mode B_A(w_x1,w_x2,w_x3,w_y)";

finish;
******************************;

start main;
** Data is collected by LF Deng 2020, with robust transformation;
use depression;
read all var _num_ into xmat;
n=nrow(xmat);
p=ncol(xmat);
ps=p*(p+1)/2;
print "N, p=" n p;
p_x=15; p_y=9;
mean_x=j(1,N,1)*xmat/N;
xmat_c=xmat-j(n,1,1)*mean_x; 
scov=xmat_c`*xmat_c/(n-1);
svar=(n-1)*vecdiag(scov)/n;
SD_in=diag(j(p,1,1)/sqrt(svar));

*starting values for ML estimation;
theta_0L={
1.006,1.057,1.202,1.085,1.110,1.193,1.288,1.335,1.071,1.221,1.399,1.340,1.234,
1.380,1.304,1.034, .916, .918, .786, .999, .947, .921, .742};
gamma_0={.066, .269, -.044};
phi_0={.451, -.182, -.187};
sig_zeta0=.202;
psi_0={
 .390, .535, .587, .379, .535, .588, .410, .612, 1.053, .823, .642, .744, .995,
 .715, .752, .174, .173, .270, .154, .305, .231, .270, .242, .310};
theta_ML=theta_0L//gamma_0//phi_0//sig_zeta0//psi_0;
q=nrow(theta_ML);
df=ps-q; print "df, q=" df q;
run DP(p, dup);
theta0=theta_ML; *starting values for ML estimation;
run ML_sem(n,p,dup,theta0,scov,T_ml,htheta_z,coef_xt,coef_yt,n_div);
if n_div=1 then do;
	print "Fisher-scoring algorithm for SEM did not converge within 300 iterations";
end;
else do;
	ncp=max(0,T_ml-df);
	RMSEA=sqrt(ncp/((n-1)*df));
	p_value=1-probchi(T_ml,df);
	print "T_ml, p_value, RMSEA=";
	print T_ml p_value RMSEA;
	id=(1:q)`;
	print "ID htheta--SE--Z";
	print ID htheta_z [format=f7.3];
	hgamma_z=htheta_z[24:26,];
	*------------------R-square under CB-SEM------------------;
	gamma=theta0[24:26];
	gamma_t=gamma`;
	phi=i(3);
	phi[2:3,1]=theta0[27:28]; phi[1,2:3]=(phi[2:3,1])`;
	phi[3,2]=theta0[29]; phi[2,3]=phi[3,2];
	*sig_zeta=theta0[30];
	*eta=gamma'\xi+zeta;
	Var_eta1=gamma_t*phi*gamma+theta0[30];
	Rsq_cb=1-theta0[30]/Var_eta1;
*---------------------------------------------------------;


	coef_x=coef_xt`;
	coef_y=coef_yt`;
	f_x=xmat_c[,1:p_x]*coef_x;
	f_y=xmat_c[,(p_x+1):p]*coef_y;
	run regress(f_x,f_y,reg_fsr, Rsq_fsr);
*----------EWC reg-------------;
	lamb_x1=theta0[1:5];
	lamb_x2=theta0[6:9];
	lamb_x3=theta0[10:15];
	lamb_y=1.0//theta0[16:23];
	c_x1=j(5,1,1)/(j(1,5,1)*lamb_x1);
	c_x2=j(4,1,1)/(j(1,4,1)*lamb_x2);
	c_x3=j(6,1,1)/(j(1,6,1)*lamb_x3);
	c_y=j(9,1,1)/(j(1,9,1)*lamb_y);
	x_c1=xmat_c[,1:5]*c_x1;
	x_c2=xmat_c[,6:9]*c_x2;
	x_c3=xmat_c[,10:15]*c_x3;
	y_c=xmat_c[,16:24]*c_y;
	x_c=x_c1||x_c2||x_c3;	
	run regress(x_c,y_c,reg_ewc, Rsq_ewc);
	print "(gamma_11//gamma_21//gamma_31)(hat\gamma, SE, z)";
	print "---CB-SEM--- ---FSR--- ---EWCR---";
	print hgamma_z[format=f7.3] reg_fsr[format=f7.3] reg_ewc[format=f7.3];
	
	xmat_c1=xmat_c[,1:5];
	xmat_c2=xmat_c[,6:9];
	xmat_c3=xmat_c[,10:15];
	ymat_c=xmat_c[,16:24];
	run PLS_A(xmat_c1,xmat_c2,xmat_c3,ymat_c, reg_plsA, reg_plsBA, Rsq_A, Rsq_BA);
	print "---plsA------plsB_A---";
	print reg_plsA[format=f7.3] reg_plsBA[format=f7.3];

*	Rsq_all=(Rsq_cb//Rsq_fsr//Rsq_ewc//Rsq_A//Rsq_BA)`;
*	print "R-square(CB-SEM,FSR,EWCR,PLSA,PLSB_A)";
*	print Rsq_all[format=f6.3];

end;

finish;
run main;