/*************************************************/
/* Written by:                                   */
/* 	                                             */
/*	G. John Geldhof                              */
/*	Graduate Student, Department of Psychology   */
/*	University of Kansas                         */
/*                                               */
/*                   AND                         */
/*                                               */
/*	Ihno Lee                                     */
/*	Graduate Student, Department of Psychology   */
/*	University of Kansas                         */
/*                                               */
/* For:                                          */
/*                                               */
/*	 HTTP://QUANT.KU.EDU                         */
/*                                               */
/*************************************************/

/*******************************************************************************************/
/*******************************************************************************************/

/*This program imports LISREL start value matrices from Y-side structural models and expands 
	them so that they can be used as starting values for a model with phantom constructs.
	The program reads in LAMBDA-Y, PSI, and ALPHA matrices, from a base model and outputs
	LAMBDA-Y, PSI, ALPHA, and BETA matrices containing the same information, but for a model
	with phantom constructs. The latent variances, for instance, are fixed to 0.0 for lower 
	order constructs and fixed to 1.0 for higher order constructs. A beta matrix is created
	that represents the latent variances. The regression path connecting the original and 
	the phantom constructs is started at a value equal to the square root of the original
	latent variance (i.e., the standard deviation of the construct). This assumes a one-to-one
	association between lower order and phantom constructs. For example a model with initial 
	constructs A,B, and C will have six constructs when phantom constructs are included: A, B,
	C, Phantom-A, Phantom-B, and Phantom-C (in that order).

	To use this program, users must enter locations for the input and output matrices, 
	and provide a few pieces of descriptive information*/

/*******************************************************************************************/
/*******************************************************************************************/


/*Indicate the location of the input matrices; only edit information in the quotes*/
         /*NOTE: TE AND TY DO NOT CHANGE DIMENSIONS AND ARE NOT REQUIRED*/

filename ly "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\strong.ly";
/*Lambda Matrix*/

filename ps "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\strong.ps";
/*Psi Matrix*/

filename al "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\strong.al";
/*Alpha Matrix*/


/*Indicate the locations and filenames for your output (this will create a set of new files)*/
filename lyout "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\phantom.ly";
/*Lambda Matrix*/

filename psout "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\phantom.ps";
/*Psi Matrix*/

filename alout "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\phantom.al";
/*Alpha Matrix*/

filename beout "C:\Documents and Settings\geldhof\Desktop\Projects\Phantom Matrices\phantom.be";
/*Beta Matrix*/



/*Tell the program how many constructs (y-side) existed when the input matrices were created.
  To do this, change the value "p2" below to read "p#", where # is the number of constructs
  present in one group of the model that created the input matrix. */

data usely; infile ly;
	input (p1-p2) (13.) @@;
	run;

/*Again, tell SAS how many latent (y-side) constructs you have by changing "p2" to "p#", 
  where # is the number of latent constructs*/

data useal; infile al;
	input (p1-p2) (13.) @@;
	run;

/*You will also need to specify how many elements there are in your psi matrix. To find out,
  specify the number of constructs in the following code and run it. This code can be run 
  by highlighting it and hitting the 'F3' key on your keyboard.*/

/*************************/
/*START # CONSTRUCTS CODE*/
/*************************/

data n_psi;
	input Number_Constructs;
	cards;
	2
	;/*Specify the number of constructs on the line above. Leaving the value at '2' specifies
	   that the initial model had two constructs*/
	run;
/*The next piece of code calculates the number of unique elements in your psi matirx*/
data n_psi; set n_psi;
	Elements = ((number_constructs*(number_constructs-1))/2)+number_constructs;
	run; 
proc print; run; /*The PROC PRINT statement prints out the number calculated above*/ 

/***********************/
/*END # CONSTRUCTS CODE*/
/***********************/

/*Now specify the number of elements in the input psi matrix, similar to what you did before*/
data useps; infile ps;
	input (p1-p3) (13.) @@;
run;

proc IML;
	/*Specify the number of constructs you want to add*/
	newconst = 2;
	/*And specify the number of groups in your model*/
	groups = 2;


/*End User Input*/


/**************************************
***************************************
********** DO NOT EDIT BELOW **********
**********    THIS POINT     **********
***************************************
**************************************/

	
	use usely; /*input the original matrices into SAS PROC IML*/
	read all into matly;
	use useal;
	read all into matal;
	use useps;
	read all into vectps;

/*Edit PS*/
zeroel = j(1,1,0); /*creates a scalar with value 0*/
do grpit = 1 to groups; /*Repeats for each group*/
	matps = j(ncol(matly),ncol(matly),0); /*creates a null psi matrix; 
											psi is read in as a vector*/
	el=1; /*the next two DO loops convert the psi vector into a symmetric matrix*/
	do i=1 to ncol(matly);
		do j=1 to i;
			matps[i,j] = vectps[grpit,el];
			matps[j,i] = vectps[grpit,el];
			el=el+1;
		end;
	end;
	dpsi = diag(matps); /*Creates a matrix containing only the diagonal
						  of the original psi matrix for standardization*/ 
	stdvect = vecdiag(matps)##.5; /*To be used when creating the beta matrix*/
	do i = 1 to newconst; /*Also used for creating the beta matrix*/
		stdvect = zeroel // stdvect;
	end;
	do i = 1 to ncol(matly); /*Takes diagonal elements of psi to the -.5 power for standardizing*/
		do j = 1 to ncol(matly);
			if dpsi[i,j] > 0 then dpsi[i,j]= dpsi[i,j] ## -.5;
		end;
	end;
	corr = dpsi * matps * dpsi; /*standardizes psi*/
	psicol0 = j(nrow(corr),1,0); 
	do i = 1 to newconst; /*Add columns of 0's to the correlation matrix for the new constructs*/
		corr = psicol0 || corr;
	end;
	psirow0 = j(1,ncol(corr),0);
	do i = 1 to newconst; /*Add rows of 0's to the correlation matrix for the new constructs*/
		corr = psirow0 // corr;
	end;
	unqcorr = ((ncol(corr) * (ncol(corr) - 1))/2) + ncol(corr); /*computes number of  elements 
				in the lower triangle of the correlation matrix*/
	corvect = j(1,unqcorr,0); /*For outputting correlations as a vector*/
	el=1;
	do i=1 to ncol(corr); /*Creates a vector of the correlation matrix*/ 
		do j=1 to i; 
			corvect[1,el] = corr[i,j];
			el=el+1;
		end;
	end;

/*Beta Matrix*/	
	bevect = j(ncol(corr), ncol(corr),0);
	do beit = 1 to newconst; /*This loopc reates an off-diagonal of the beta matrix that
								corresponds to the original latent standard deviations*/
		bevect[beit,newconst + beit] = stdvect[newconst + beit,1];
	end;
	if grpit = 1 then outpsi = corvect; /*Outputs and concactenates across groups*/
	if grpit = 1 then outbe = bevect;
	if grpit >1 then outpsi = outpsi // corvect;
	if grpit >1 then outbe = outbe // bevect;

end;

/*Edit LY*/
	vect0 = j(nrow(matly),1,0);
	outly = matly;
	do i = 1 to newconst;
		outly = outly || vect0; /*Adds columns of 0's to LY to accomodate new 
								  latent constructs*/
	end;
	
/*Edit AL*/
	single0 = j(nrow(matal),1,0);
	outal = matal;
	do i = 1 to newconst;
		outal = outal || single0; /*Adds columns of 0's to AL to accomodate new 
								    latent constructs*/
	end;

	print outpsi; /*prints edited matrices*/
	print outbe;
	print outly;
	print outal;

/*Output*/
	create newps FROM outpsi; /*Creates a new SAS data set called 'newps'*/
	append FROM outpsi;	/*Outputs the data in matrix 'outpsi' into the newly created file*/

	create newbe FROM outbe;
	append FROM outbe;

	create newly FROM outly;
	append FROM outly;

	create newal FROM outal;
	append FROM outal;
		
quit;
/*Outputting the new SAS datasets into external files*/
data newps; set newps; /*Adds an additional (blank) variable to the data set. this is neccessary
						 for outputting the data to an external file*/
	blank = ' ';
run;
data newbe; set newbe;
	blank = ' ';
run;
data newly; set newly;
	blank = ' ';
run;
data newal; set newal;
	blank = ' ';
run;
data test; set newps; file psout;
	put col1--blank; /*Outputs all data between the first column and the blank variable 
					   into the new file*/
run;
data test; set newbe; file beout;
	put col1--blank;
run;
data test; set newly; file lyout;
	put col1--blank;
run;
data test; set newal; file alout;
	put col1--blank;
run;