**********************************************************************************
**This is to generate 500 datasets with continuous distribution         **
**under the condition of:                                                                        **
**sample size =200,                                                                              **
**loading=0.8 for all items,                                                                    **
**underlying distribution=normal,                                                          **
**number of scale points=2,                                                                 **
**and the threshold for creating the categories=0                              **  
*********************************************************************************;

*********************************************************************************
**Note that the files "mlm.inp", "wlsmv.inp", and "replist.dat"             **
**containing the name of 500 datasets (rep1.dat-rep500.dat)          **
**needed to be created before running the codes.                           **
*********************************************************************************;

libname SEM 'd:\GY'; 
%let reps=500; 
%Macro GY;
%do t=1 %to &reps; 

proc iml;
 CORR={1 .64 .64 .64 .64 .64 .64 .64,
       .64 1 .64 .64 .64 .64 .64 .64,
       .64 .64 1 .64 .64 .64 .64 .64,
       .64 .64 .64 1 .64 .64 .64 .64,
       .64 .64 .64 .64 1 .64 .64 .64,
       .64 .64 .64 .64 .64 1 .64 .64,
       .64 .64 .64 .64 .64 .64 1 .64,
       .64 .64 .64 .64 .64 .64 .64 1};

X=rannor(J(200,8,0));
root=root(CORR);
X=X*root;  

x1=X[,1];X2=X[,2];X3=X[,3];x4=X[,4];
x5=X[,5];x6=X[,6];x7=X[,7];x8=X[,8];

data=X1||X2||x3||X4||x5||X6||x7||X8; 

create dataset from data[colname={x1 x2 x3 x4 x5 x6 x7 x8}];
append from data; 

data temp;set dataset;
	 if x1<=0 then y1=0; else y1=1;
     if x2<=0 then y2=0; else y2=1;
     if x3<=0 then y3=0; else y3=1;
     if x4<=0 then y4=0; else y4=1;
     if x5<=0 then y5=0; else y5=1;
     if x6<=0 then y6=0; else y6=1;
     if x7<=0 then y7=0; else y7=1;
     if x8<=0 then y8=0; else y8=1;	 
run;

data _null_;
set temp;
file "d:\GY\rep&t..dat";
put y1 y2 y3 y4 y5 y6 y7 y8;
run;

%end;
%Mend GY;
%GY

options noxwait xsync;
x 'cd c:\' ;
x 'cd C:\Program Files\Mplus' ;
x "mplus D:\GY\mlm.inp D:\GY\mlm.out" ;

options noxwait xsync;
x 'cd c:\' ;
x 'cd C:\Program Files\Mplus' ;
x "mplus D:\GY\wlsmv.inp D:\GY\wlsmv.out" ;

data temp; 
infile 'D:\GY\wlsmv.dat';
input 
rep tau1 tau2 tau3 tau4 tau5 tau6 tau7 tau8
y1 y2 y3 y4 y5 y6 y7 y8
sdtau1-sdtau8 sdy1-sdy8
chisq df p_chisq CFI TLI npara RMSEA RMSEALO RMSEAHI RMSEAP WRMR;
run;

%MACRO nonlinear;
%DO r = 1 %to 500;
data input;
set temp  (firstobs=&r obs=&r);

data input1;set input (keep=tau1 tau2 tau3 tau4 tau5 tau6 tau7 tau8);
run;
data input2;set input (keep=y1 y2 y3 y4 y5 y6 y7 y8);
run;

proc iml;
use input1; 
read all var _num_ into TTHRESH;
run;
use input2; 
read all var _num_ into TLOAD;
run;

THRESH=T(TTHRESH);
LOAD=T(TLOAD);

FACCOR={1};
NTHRESH=Ncol(thresh);
NCAT=NTHRESH+1;
NITEM=Nrow(LOAD);
NFACT=Ncol(LOAD);
POLYR=LOAD*FACCOR*T(LOAD);
do j=1 to NITEM;
POLYR[j,j]=1;
end;
sumnum=0;
addden=0;
do j=1 to NITEM;
do jp=1 to NITEM;
sumprobn2=0;
addprobn2=0;
do c=1 to NTHRESH;
do cp=1 to NTHRESH;
sumrvstar=0;
do k=1 to NFACT;
do kp=1 to NFACT;
sumrvstar=sumrvstar+LOAD[j,k]*LOAD[jp,kp]*FACCOR[k,kp];
end;
end;
sumprobn2=sumprobn2+probbnrm(THRESH[j,c],THRESH[jp,cp],sumrvstar);
addprobn2=addprobn2+probbnrm(THRESH[j,c],THRESH[jp,cp],POLYR[j,jp]);
end;
end;
sumprobn1=0;
sumprobn1p=0;
do cc=1 to NTHRESH;
sumprobn1=sumprobn1+CDF('NORMAL',THRESH[j,cc]);
sumprobn1p=sumprobn1p+CDF('NORMAL',THRESH[jp,cc]);
end;
sumnum=sumnum+(sumprobn2-sumprobn1*sumprobn1p);
addden=addden+(addprobn2-sumprobn1*sumprobn1p);
end;
end;
reliab=sumnum/addden;

data=&r||sumnum||addden||reliab; 

create nonlinear from data[colname={rep num den reliab}];
append from data; 

proc append base=sem.wlsmv_reliab data=nonlinear force;
run;

%end;
%mend nonlinear;
%nonlinear;