/*
******************************************************************************************;
*
*CTANEST1.sas
*By John R. Hipp and Daniel J. Bauer;
*
*Updated June 5, 2003
*
*Adding randomization of tetrads (does this work for nested models???);
*
*This compares nested models using the Tetrad test;
*
*If using categorical data, use CTAFILE.mac to prepare the ACM file;
*
******************************************************************************************;
*/


%macro ctanest1(N=, vars=, nesttest=0, pchor=0, lisrel=0, mplus=0, SAMPMAT1=, IMPMAT1B=, 
				IMPMAT2B=, lowdiag=0, reps=1);

%if &lowdiag = 1 %then %do;
	%let fullmat = %eval((&vars + 1)*(&vars)/2);
%end;
%else %if &lowdiag = 0 %then %do;
	%let fullmat = %eval(&vars*&vars);
%end;

%if &pchor = 0 %then %do;
data t;
	array t(&fullmat) (&SAMPMAT1);
run;
%end;

data imp1;
	array imp1(&fullmat) (&IMPMAT1b);
run;

/* If performing nested tests, read in the second model imp cov, otherwise just use the same one */
%if &nesttest = 1 %then %do;
	data imp2;
		array imp2(&fullmat) (&IMPMAT2b);
	run;
%end;
%else %do;
	data imp2;
		array imp2(&fullmat) (&IMPMAT1b);
	run;
%end;



proc iml;


START TETRAD;

	IF &nesttest = 0 THEN DO;
		tetlabel = initlabel_nonnest;
		tetres = initres_nonnest;
	END;
	ELSE DO;
		tetlabel = initlabel;
		tetres = initres;
	END;

	/* indexing tetrads; */
	id = 1;
	DO i = 2 to NROW(tetres);
		id = id//i;
	END;

	**Column 9 is id;
	tetlabel = tetlabel||id;

	IF &reps > 1 THEN DO;
		*RANDOMIZE the tetrads;
		ran = J(NROW(tetlabel),1,0);
		DO i = 1 to NROW(tetlabel);
			ran[i,1]=uniform(0);
		END;
		tetlabel = tetlabel||ran;
		tetlabel2 = tetlabel;
		tetlabel[rank(tetlabel[,10]),]=tetlabel2;
		tetres = tetres||ran;
		tetres2 = tetres;
		tetres[rank(tetres[,2]),]=tetres2;
		/* Now drop the random column from tetlabel */
		tetlabel = tetlabel[,1:9];
	END;



/*
***********************************************************************;
*Now loop through and find the vanishing tetrads for the model with fewer implied;
***********************************************************************;
*Loop here only if there are tetrads; */
IF &nesttest = 1 & zerotet = 0 THEN DO;

	/* The tenth column contains the residuals */
	tetlabel = tetlabel||tetres[,1];

	/* find unique covariances among implied tetrads */
	/* Note that the order of the tetrads gets flipped here to match with the order that 
	LISREL and M-Plus give in their ACM's 
	(do this by placing the j+1 element in front of the j element); */
	tetlabel2 = J(NROW(tetlabel),4,'aaaa');
	DO i = 1 to NROW(tetlabel);
		k = 0;
		DO j = 1 to 7 BY 2;
			k = k+1;
			IF tetlabel[i,j] < 10 & tetlabel[i,j+1] < 10 THEN 
				tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			ELSE DO;
				IF tetlabel[i,j] >= 10 & tetlabel[i,j+1]>=10 THEN 
					tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
				ELSE IF tetlabel[i,j] >= 10 THEN
						tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
					ELSE tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			END;
		END;
	END;
	tetcov = UNIQUE(SHAPE(tetlabel2,1)); /* row vector of unique covariances; */

	/*
	**************************************************************************;
    	  construct derivative matrix: d 
	**************************************************************************;
	*/
	D = J(NCOL(tetcov),NROW(tetlabel2),0);
	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NROW(tetlabel2);
			IF tetcov[,i] = tetlabel2[j,1]
				THEN D[i,j] = IMPMAT2[tetlabel[j,3],tetlabel[j,4]];
			IF tetcov[,i] = tetlabel2[j,2]
				THEN D[i,j] = IMPMAT2[tetlabel[j,1],tetlabel[j,2]];
			IF tetcov[,i] = tetlabel2[j,3]
				THEN D[i,j] = -1*IMPMAT2[tetlabel[j,7],tetlabel[j,8]];
			IF tetcov[,i] = tetlabel2[j,4]
				THEN D[i,j] = -1*IMPMAT2[tetlabel[j,5],tetlabel[j,6]];
		END;
	END;

	/* compute asymptotic covariance matrix of the IMPLIED covariances */
	SigmaS = J(NCOL(tetcov),NCOL(tetcov),0);
	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NCOL(tetcov);
			e = NUM(SUBSTR(tetcov[,i],1,2));
			f = NUM(SUBSTR(tetcov[,i],3,2));
			g = NUM(SUBSTR(tetcov[,j],1,2));
			h = NUM(SUBSTR(tetcov[,j],3,2));
			/* This calculates it for a correlation matrix; */
			IF corr = 1 THEN DO;
				SigmaS[i,j] = (1/2)
                	   	* (IMPMAT2[e,f]*IMPMAT2[g,h])
                   		* (IMPMAT2[e,g]**2 + IMPMAT2[e,h]**2 + IMPMAT2[f,g]**2 + IMPMAT2[f,h]**2)
                   		+ IMPMAT2[e,g]*IMPMAT2[f,h] + IMPMAT2[e,h]*IMPMAT2[f,g]
                   		- IMPMAT2[e,f]*(IMPMAT2[f,g]*IMPMAT2[f,h]+IMPMAT2[e,g]*IMPMAT2[e,h])
                   		- IMPMAT2[g,h]*(IMPMAT2[f,g]*IMPMAT2[e,g]+IMPMAT2[f,h]*IMPMAT2[e,h]);
			END;
			/* This calculates it for a covariance matrix; */
			ELSE DO;
				SigmaS[i,j] = IMPMAT2[e,g]*IMPMAT2[f,h] + IMPMAT2[e,h]*IMPMAT2[f,g];
			END;
		END;
	END;

	SigmaTemp = D`*SigmaS*D; 

	/* Add an 11th column for marking tetrads to be dropped; */
	tetlabel = tetlabel||J(NROW(tetlabel),1,0);

	/*
	*********************************************************************;
	*Using the SWEEP operator;
	*********************************************************************;
	*/
	/* ***Don't go into this loop if only one vanishing tetrad; */
	IF nrow(tetlabel) > 1 THEN DO;

		temp33 = sweep(SigmaTemp);
		temp33 = temp33||J(NROW(temp33),1,0);
 
		DO i = 1 to NROW(temp33);
			IF sum(abs(temp33[i,])) <> 0 THEN tetlabel[i,11]=1;
		END;

		/* Now drop the zero rows and columns (singular); */
		temp5 = J(1,NCOL(temp33),0);

		DO i = 1 to NROW(temp33);
			IF sum(abs(temp33[i,])) <> 0 THEN DO;
				temp5 = temp5 // temp33[i,1:NCOL(temp33)];
			END;
		END;
		temp2 = J(NROW(temp5),1,0);
		DO i = 1 to NCOL(temp5);
			IF sum(abs(temp5[,i])) <> 0 THEN DO;
				temp2 = temp2 || temp5[1:NROW(temp5),i];
			END;
		END;

		SigmaTemp2 = temp2[2:NROW(temp2),2:NCOL(temp2)];

		temp7 = J(1,11,0);
		temp8 = J(1,11,0);

		/* Now find the redundant tetrads; */
		DO i = 1 to NROW(tetlabel);
			IF tetlabel[i,11]=1 then DO;
				temptet= tetlabel[i,];
				temp7 = temp7//temptet;
			END;
			IF tetlabel[i,11]=0 then DO;
				temptet2= tetlabel[i,];
				temp8 = temp8//temptet2;
			END;
		END;
		IF NROW(temp7) > 1 THEN DO;
			temp7a = temp7[2:NROW(temp7),1:8];
			temp7b = temp7[2:NROW(temp7),10];
		END;
		IF NROW(temp8) > 1 THEN DO;
			temp8a = temp8[2:NROW(temp8),1:8];
			temp8b = temp8[2:NROW(temp8),10];
		END;

		/* These are the tetrads in the smaller model, and their residuals; */
		tetlabel_small = temp7a;
		tetres_small = temp7b;

		/*
		*******************************************************************************;
		*END of loop--now reorder the tetrads so that we're pulling the same ones in the larger model;
		*******************************************************************************;
		*/
		tetlabel = J(NROW(tetlabel),8,0);

		/* Order is nonredundant---not in model---redundant; */
		tetlabel = temp7a//temp3//temp8a;
		tetres = temp7b//temp4//temp8b;
	END;
	ELSE DO;
		tetlabel_small = tetlabel[1,1:8];
		tetres_small = tetlabel[1,10];
		tetlabel = tetlabel_small;
		tetres = tetres_small;
	END;
END;
FREE Temp33 Temp5 Temp2 Temp7 Temp8 Sigmatemp2 Temptet Temptet2 Temp7a Temp7b Temp8a Temp8b;
/*
***********************************************************************;
*Now loop through a SECOND time and find the vanishing tetrads for the model with more implied;
***********************************************************************;
*/
/* This used to select on some instances---doesn't appear necessary anymore; */

	IF &nesttest = 1 THEN DO;
		/* indexing tetrads; */
		id = 1;
		DO i = 2 to NROW(tetres);
			id = id//i;
		END;
		/* Column 9 is id; */
		tetlabel = tetlabel||id;
	END;
 


	/* find unique covariances among implied tetrads */
	tetlabel2 = J(NROW(tetlabel),4,'aaaa');
	DO i = 1 to NROW(tetlabel);
		k = 0;
		DO j = 1 to 7 BY 2;
			k = k+1;
			IF tetlabel[i,j] < 10 & tetlabel[i,j+1] < 10 THEN 
				tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			ELSE DO;
				IF tetlabel[i,j] >= 10 & tetlabel[i,j+1]>=10 THEN 
					tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
				ELSE IF tetlabel[i,j] >= 10 THEN
						tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
					ELSE tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			END;
		END;
	END;

	tetcov = UNIQUE(SHAPE(tetlabel2,1)); /* row vector of unique covariances; */

	/*
	**************************************************************************;
    	  construct derivative matrix: d 
	**************************************************************************;
	*/
	D = J(NCOL(tetcov),NROW(tetlabel2),0);

	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NROW(tetlabel2);
			IF tetcov[,i] = tetlabel2[j,1]
				THEN D[i,j] = IMPMAT[tetlabel[j,3],tetlabel[j,4]];
			IF tetcov[,i] = tetlabel2[j,2]
				THEN D[i,j] = IMPMAT[tetlabel[j,1],tetlabel[j,2]];
			IF tetcov[,i] = tetlabel2[j,3]
				THEN D[i,j] = -1*IMPMAT[tetlabel[j,7],tetlabel[j,8]];
			IF tetcov[,i] = tetlabel2[j,4]
				THEN D[i,j] = -1*IMPMAT[tetlabel[j,5],tetlabel[j,6]];
		END;
	END;

	/* compute asymptotic covariance matrix of the IMPLIED covariances */
	SigmaS = J(NCOL(tetcov),NCOL(tetcov),0);
	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NCOL(tetcov);
			e = NUM(SUBSTR(tetcov[,i],1,2));
			f = NUM(SUBSTR(tetcov[,i],3,2));
			g = NUM(SUBSTR(tetcov[,j],1,2));
			h = NUM(SUBSTR(tetcov[,j],3,2));
			/* This calculates it for a correlation matrix; */
			IF corr = 1 THEN DO;
			SigmaS[i,j] = (1/2)
            	    * (IMPMAT[e,f]*IMPMAT[g,h])
                   	* (IMPMAT[e,g]**2 + IMPMAT[e,h]**2 + IMPMAT[f,g]**2 + IMPMAT[f,h]**2)
                   	+ IMPMAT[e,g]*IMPMAT[f,h] + IMPMAT[e,h]*IMPMAT[f,g]
                   	- IMPMAT[e,f]*(IMPMAT[f,g]*IMPMAT[f,h]+IMPMAT[e,g]*IMPMAT[e,h])
                   	- IMPMAT[g,h]*(IMPMAT[f,g]*IMPMAT[e,g]+IMPMAT[f,h]*IMPMAT[e,h]);
			END;
			/* This calculates it for a covariance matrix; */
			ELSE DO;
				SigmaS[i,j] = IMPMAT[e,g]*IMPMAT[f,h] + IMPMAT[e,h]*IMPMAT[f,g];
			END;
		END;
	END;

	tetlabel = tetlabel||J(NROW(tetlabel),1,0);
	SigmaTemp = D`*SigmaS*D; 
	FREE SigmaS D;
	/*
	*********************************************************************;
	*Using the SWEEP operator;
	*********************************************************************;
	*/
	/* ***Don't go into this loop if only one vanishing tetrad; */
	IF nrow(tetlabel) > 1 THEN DO;

		temp33 = sweep(SigmaTemp);
		temp33 = temp33||J(NROW(temp33),1,0);

		DO i = 1 to NROW(temp33);
			IF sum(abs(temp33[i,])) <> 0 THEN tetlabel[i,10]=1;
		END;

		/* Now drop the zero rows and columns (singular); */
		temp5 = J(1,NCOL(temp33),0);

		DO i = 1 to NROW(temp33);
			IF sum(abs(temp33[i,])) <> 0 THEN DO;
				temp5 = temp5 // temp33[i,1:NCOL(temp33)];
			END;
		END;
		temp2 = J(NROW(temp5),1,0);
		DO i = 1 to NCOL(temp5);
			IF sum(abs(temp5[,i])) <> 0 THEN DO;
				temp2 = temp2 || temp5[1:NROW(temp5),i];
			END;
		END;

		SigmaTemp2 = temp2[2:NROW(temp2),2:NCOL(temp2)];

		/* Keep the nonredundant tetrads; */
		DO i = 1 to NROW(tetlabel);
			IF tetlabel[i,10]=1 then DO;
				temptet= tetlabel[i,];
				tetlab4 = tetlab4//temptet;
				tetlab5 = tetlab5//tetres[i,1];
			END;
		END;
		tetlabel = tetlab4[,1:9];
		tetres = tetlab5;
	END;
	ELSE DO;
		tetlabel = tetlabel[,1:9];
	END;
	FREE Sigmatemp Temp33 Temp5 Temp2 SigmaTemp2 Temptet Tetlab4 Tetlab5;

	/*
	*********************************************************************;
	*Now use the SWEEP results on the original covariance matrix of tetrads
	*to get the model fit for the larger model;
	*********************************************************************;
	*/
	/* find unique covariances among non-redundant tetrads */
	tetlabel2 = J(NROW(tetlabel),.5*(NCOL(tetlabel)-1),'aaaa');
	DO i = 1 to NROW(tetlabel);
		k = 0;
		DO j = 1 to 7 BY 2;
			k = k+1;
			IF tetlabel[i,j] < 10 & tetlabel[i,j+1] < 10 THEN 
				tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			ELSE DO;
				IF tetlabel[i,j] >= 10 & tetlabel[i,j+1]>=10 THEN 
					tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
				ELSE IF tetlabel[i,j] >= 10 THEN
						tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
					ELSE tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			END;
		END;
	END;

	tetcov = UNIQUE(SHAPE(tetlabel2,1)); /* row vector of unique covariances; */

	/*
	**************************************************************************;
     	 construct derivative matrix: d 
	**************************************************************************;
	*/
	D = J(NCOL(tetcov),NROW(tetlabel2),0);
	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NROW(tetlabel2);
			IF tetcov[,i] = tetlabel2[j,1]
				THEN D[i,j] = SAMPMAT[tetlabel[j,3],tetlabel[j,4]];
			IF tetcov[,i] = tetlabel2[j,2]
				THEN D[i,j] = SAMPMAT[tetlabel[j,1],tetlabel[j,2]];
			IF tetcov[,i] = tetlabel2[j,3]
				THEN D[i,j] = -1*SAMPMAT[tetlabel[j,7],tetlabel[j,8]];
			IF tetcov[,i] = tetlabel2[j,4]
				THEN D[i,j] = -1*SAMPMAT[tetlabel[j,5],tetlabel[j,6]];
		END;
	END;

	/* compute asymptotic covariance matrix of the covariances -- normal theory estimator;
	Skip this if bringing in an outside polychoric matrix; */

	IF &lisrel = 0 & &mplus = 0 THEN DO;
		SigmaS = J(NCOL(tetcov),NCOL(tetcov),0);
		DO i = 1 to NCOL(tetcov);
			DO j = 1 to NCOL(tetcov);
				e = NUM(SUBSTR(tetcov[,i],1,2));
				f = NUM(SUBSTR(tetcov[,i],3,2));
				g = NUM(SUBSTR(tetcov[,j],1,2));
				h = NUM(SUBSTR(tetcov[,j],3,2));
				/* This calculates it for a correlation matrix; */
				IF corr = 1 THEN DO;
				SigmaS[i,j] = (1/2)
                    	* (SAMPMAT[e,f]*SAMPMAT[g,h])
                    	* (SAMPMAT[e,g]**2 + SAMPMAT[e,h]**2 + SAMPMAT[f,g]**2 + SAMPMAT[f,h]**2)
                    	+ SAMPMAT[e,g]*SAMPMAT[f,h] + SAMPMAT[e,h]*SAMPMAT[f,g]
                    	- SAMPMAT[e,f]*(SAMPMAT[f,g]*SAMPMAT[f,h]+SAMPMAT[e,g]*SAMPMAT[e,h])
                    	- SAMPMAT[g,h]*(SAMPMAT[f,g]*SAMPMAT[e,g]+SAMPMAT[f,h]*SAMPMAT[e,h]);
				END;
				/* This calculates it for a covariance matrix; */
				ELSE DO;
				SigmaS[i,j] = SAMPMAT[e,g]*SAMPMAT[f,h] + SAMPMAT[e,h]*SAMPMAT[f,g];
				END;
			END;
		END;
	END;
	ELSE DO;
		/* Read in the ACM from LISREL or M-Plus; */
		temp1 = J(1,NCOL(acmcov),0);
		DO i = 1 TO NCOL(acmcov);
			DO j = 1 TO NCOL(tetcov);
				IF acmcov[1,i] = tetcov[1,j] THEN temp1[1,i]=i;
			END;
		END;
		acmdrop = J(1,1,0);
		DO i = 1 TO NCOL(acmcov);
			IF temp1[1,i] > 0 THEN acmdrop = acmdrop||temp1[1,i];
		END;
		acmdrop2 = acmdrop[,2:NCOL(acmdrop)];

		/* Use the ACM brought in at the top; */
		SigmaS = symmat[acmdrop2,acmdrop2];
		FREE temp1 acmdrop acmdrop2;
	END;


	/* Perform Test */
	SigmaT = D`*SigmaS*D; 
	df = NROW(SigmaT);
	Chi = &N*(tetres`*INV(SigmaT)*tetres);
	p = 1-PROBCHI(Chi,df);
		IF &nesttest = 1 THEN PRINT "Tetrad Test for Model with more vanishing tetrads";
		ELSE PRINT "Tetrad Test for the Model";
	PRINT Chi[label='Chi-Square'] df[label='df'] p[label='p-value']; 

	Chibig = chi;
	dfbig = df;
	pbig = p;

	/* This collects the randomization results by replication if testing just one model */
	IF &reps > 1 & &nesttest = 0 THEN DO;
		IF r = 1 THEN DO;
			results = r||chibig||dfbig||pbig; 
			clabel = {'replication' '      chi' '  df' '  p'};
		END;
		ELSE DO;
			represult = r||chibig||dfbig||pbig;
			results = results//represult;
		END;
	END;

	FREE SigmaT SigmaS D;

/*
*****************************************************************;
*Now run the smaller model;
*****************************************************************;
*/
IF &nesttest = 1 & zerotet = 0 THEN DO;

	tetlabel = tetlabel_small;
	tetres = tetres_small;

	/* find unique covariances among non-redundant tetrads */
	tetlabel2 = J(NROW(tetlabel),.5*(NCOL(tetlabel)),'aaaa');
	DO i = 1 to NROW(tetlabel);
		k = 0;
		DO j = 1 to 7 BY 2;
			k = k+1;
			IF tetlabel[i,j] < 10 & tetlabel[i,j+1] < 10 THEN 
				tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			ELSE DO;
				IF tetlabel[i,j] >= 10 & tetlabel[i,j+1]>=10 THEN 
					tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
				ELSE IF tetlabel[i,j] >= 10 THEN
						tetlabel2[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
					ELSE tetlabel2[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			END;
		END;
	END;
	tetcov = UNIQUE(SHAPE(tetlabel2,1)); /* row vector of unique covariances; */

	/*
	**************************************************************************;
    	  construct derivative matrix: d 
	**************************************************************************;
	*/
	D = J(NCOL(tetcov),NROW(tetlabel2),0);
	DO i = 1 to NCOL(tetcov);
		DO j = 1 to NROW(tetlabel2);
			IF tetcov[,i] = tetlabel2[j,1]
				THEN D[i,j] = SAMPMAT[tetlabel[j,3],tetlabel[j,4]];
			IF tetcov[,i] = tetlabel2[j,2]
				THEN D[i,j] = SAMPMAT[tetlabel[j,1],tetlabel[j,2]];
			IF tetcov[,i] = tetlabel2[j,3]
				THEN D[i,j] = -1*SAMPMAT[tetlabel[j,7],tetlabel[j,8]];
			IF tetcov[,i] = tetlabel2[j,4]
				THEN D[i,j] = -1*SAMPMAT[tetlabel[j,5],tetlabel[j,6]];
		END;
	END;


	/* compute asymptotic covariance matrix of the covariances -- normal theory estimator;
	Skip this if bringing in an outside polychoric matrix; */

	IF &lisrel = 0 & &mplus = 0 THEN DO;
		SigmaS = J(NCOL(tetcov),NCOL(tetcov),0);
		DO i = 1 to NCOL(tetcov);
			DO j = 1 to NCOL(tetcov);
				e = NUM(SUBSTR(tetcov[,i],1,2));
				f = NUM(SUBSTR(tetcov[,i],3,2));
				g = NUM(SUBSTR(tetcov[,j],1,2));
				h = NUM(SUBSTR(tetcov[,j],3,2));
				/* This calculates it for a correlation matrix; */
				IF corr = 1 THEN DO;
				SigmaS[i,j] = (1/2)
                    	* (SAMPMAT[e,f]*SAMPMAT[g,h])
                    	* (SAMPMAT[e,g]**2 + SAMPMAT[e,h]**2 + SAMPMAT[f,g]**2 + SAMPMAT[f,h]**2)
                    	+ SAMPMAT[e,g]*SAMPMAT[f,h] + SAMPMAT[e,h]*SAMPMAT[f,g]
                    	- SAMPMAT[e,f]*(SAMPMAT[f,g]*SAMPMAT[f,h]+SAMPMAT[e,g]*SAMPMAT[e,h])
                    	- SAMPMAT[g,h]*(SAMPMAT[f,g]*SAMPMAT[e,g]+SAMPMAT[f,h]*SAMPMAT[e,h]);
				END;
				/* This calculates it for a covariance matrix; */
				ELSE DO;
				SigmaS[i,j] = SAMPMAT[e,g]*SAMPMAT[f,h] + SAMPMAT[e,h]*SAMPMAT[f,g];
				END;
			END;
		END;
	END;
	ELSE DO;
		/* Read in the ACM from LISREL or M-Plus; */
		temp1 = J(1,NCOL(acmcov),0);
		DO i = 1 TO NCOL(acmcov);
			DO j = 1 TO NCOL(tetcov);
				IF acmcov[1,i] = tetcov[1,j] THEN temp1[1,i]=i;
			END;
		END;
		acmdrop = J(1,1,0);
		DO i = 1 TO NCOL(acmcov);
			IF temp1[1,i] > 0 THEN acmdrop = acmdrop||temp1[1,i];
		END;
		acmdrop2 = acmdrop[,2:NCOL(acmdrop)];

		/* Use the ACM brought in at the top; */
		SigmaS = symmat[acmdrop2,acmdrop2];
	END;

	/* Perform Test */
	SigmaT = D`*SigmaS*D; 
	df = NROW(SigmaT);
	Chi = &N*(tetres`*INV(SigmaT)*tetres);
	p = 1-PROBCHI(Chi,df);
	PRINT "Tetrad Test for Model with fewer Vanishing tetrads";
	PRINT Chi[label='Chi-Square'] df[label='df'] p[label='p-value']; 
	Chitest = Chibig - Chi;
	dftest = dfbig - df;
	ptest = 1-PROBCHI(Chitest,dftest);
	PRINT "Nested Tetrad Test for two models";
	PRINT Chitest[label='Chi-Square'] dftest[label='df'] ptest[label='p-value']; 

	/* This collects the randomization results by replication */
	IF &reps > 1 THEN DO;
		IF r = 1 THEN DO;
			results = r||chibig||dfbig||pbig||chi||df||p||chitest||dftest||ptest; 
			clabel = {'  reps' ' big-chi' 'big-df' '  big-p' 'small-chi' 'small-df' ' small-p' 'test-chi' 'test-df' ' test-p'};
		END;
		ELSE DO;
			represult = r||chibig||dfbig||pbig||chi||df||p||chitest||dftest||ptest;
			results = results//represult;
		END;
	END;


END;

FINISH;



/********************************************************************************/
/*  MODULE:  MAIN    									                        */
/*  LAST MODIFIED:  6/5/2003                                                   */
/*  																		  	*/
/*  Gives requisite matrices, sample size, loops through Tetrad module			*/
/********************************************************************************/

START;

/* Setting a dummy to check for models implying NO vanishing tetrads; */
zerotet = 0;
N = &N;
notnest = 99;
matswitch = 99;

IF &lisrel = 1 | &mplus = 1 THEN DO;
	use acm1;
	read all into c1;
	
	s=J(1,1,0);                              /*Set up dummy row for concatenation */
	DO i=1 to ncol(c1);
		IF c1[i] =0 then
    	DO;
		
		END;                        /*Do nothing on these zero-ones, but don't advance counter */
		ELSE IF c1[i] =1 then
    	DO;

		END;
    	ELSE DO;
		s = s//c1[i];
		END;
	END;
	FREE c1;
	s = s[2:NROW(s),];				/*Removing dummy row; */

	/* Multiply this by N for the M-Plus ACM; */
	IF &mplus = 1 THEN DO;
		s2 = s*&N;
		s = s2;
	END;
		
	/*Now convert this vector into a square ACM matrix to use in CTA; */
	ssize = (&vars*(&vars-1)/2);
	acm=J(ssize,ssize,0);
	ct = 1;
	DO i=1 to nrow(acm);
    	DO j=1 to nrow(acm);
        	IF i < j then DO;
        		acm[i,j] = 0;
        	END;
    		ELSE DO;
        		acm[i,j] = s[ct];
        		ct = ct + 1;
        	END;
    	END;
	END;
	/*Symmetrize this matrix; */
 	symmat=acm+acm`;
	DO i=1 to nrow(symmat);
    	DO j=1 to nrow(symmat);
			IF i=j then symmat[i,j]= symmat[i,j]/ 2;
		END;
	END;
END;
FREE ssize acm;

IF &pchor = 1 THEN DO;
	/*Bringing in the file of the Polychoric Matrix; */
	use pchor1;
	read all into d1a;

	/*Now convert this vector into a square matrix to use in CTA; */
	pcm=J(&vars,&vars,0);
	ct = 1;
	DO i=1 to nrow(pcm);
    	DO j=1 to nrow(pcm);
        	IF i < j then DO;
        		/*Do nothing, since already a 0 value; */
        	END;
    		ELSE DO;
        		pcm[i,j] = d1a[ct];
        		ct = ct + 1;
        	END;
    	END;
	END;

	/*Symmetrize this matrix; */
 	SAMPMAT=pcm+pcm`;
	DO i=1 to nrow(SAMPMAT);
    	DO j=1 to nrow(SAMPMAT);
			IF i=j then SAMPMAT[i,j]= SAMPMAT[i,j]/ 2;
		END;
	END;
FREE d1a pcm;
END;
ELSE DO;

	/* Convert the SAMPMAT1 vector into a square matrix to use in CTA; */
	/* Symmetrize this matrix, if not already symmetric; */

	use t;
	read all into SAMPMAT1;

	IF NCOL(SAMPMAT1) = (&vars*(&vars+1)/2) THEN DO;
	/* IF lowdiag = 1 THEN DO; */
		SAMPMAT=J(&vars,&vars,0);
		ct = 1;
		DO i=1 to nrow(SAMPMAT);
    		DO j=1 to nrow(SAMPMAT);
        		IF i < j then DO;
        			SAMPMAT[i,j] = 0;
        		END;
    			ELSE DO;
        			SAMPMAT[i,j] = SAMPMAT1[ct];
        			ct = ct + 1;
        		END;
    		END;
		END;
		tempmat=sampmat+sampmat`;
		DO i=1 to nrow(tempmat);
    		DO j=1 to nrow(tempmat);
				IF i=j then tempmat[i,j]= tempmat[i,j]/ 2;
			END;
		END;
		SAMPMAT = tempmat;
		FREE tempmat;
	END;
	ELSE IF NCOL(SAMPMAT1) = &vars * &vars THEN DO;
		/* Convert this vector into a square matrix to use in CTA; */
		SAMPMAT=J(&vars,&vars,0);
		ct = 1;
		DO i=1 to nrow(SAMPMAT);
    		DO j=1 to nrow(SAMPMAT);
       			SAMPMAT[i,j] = SAMPMAT1[ct];
       			ct = ct + 1;
    		END;
		END;
	END;
	ELSE DO;
		PRINT "The sample matrix entered does not conform to the # of variables specified";
	END;
END;



/*Model implied cov matrix for the model with MORE vanishing tetrads;
Convert the IMPMAT1B vector into a square matrix to use in CTA;
Symmetrize this matrix, if not already symmetric; */

use imp1;
read all into IMPMAT1B; 

IF NCOL(IMPMAT1B) = (&vars*(&vars+1)/2) THEN DO;
/*IF lowdiag = 1 THEN DO; */
	IMPMAT1A=J(&vars,&vars,0);
	ct = 1;
	DO i=1 to nrow(IMPMAT1A);
   		DO j=1 to nrow(IMPMAT1A);
       		IF i < j then DO;
       			IMPMAT1A[i,j] = 0;
       		END;
   			ELSE DO;
       			IMPMAT1A[i,j] = IMPMAT1B[ct];
       			ct = ct + 1;
       		END;
   		END;
	END;
	tempmat=IMPMAT1A+IMPMAT1A`;
	DO i=1 to nrow(tempmat);
   		DO j=1 to nrow(tempmat);
			IF i=j then tempmat[i,j]= tempmat[i,j]/ 2;
		END;
	END;
	IMPMAT1A = tempmat;
END;
ELSE IF NCOL(IMPMAT1B) = &vars * &vars THEN DO;
	/*Convert this vector into a square matrix to use in CTA; */
	IMPMAT1A=J(&vars,&vars,0);
	ct = 1;
	DO i=1 to nrow(IMPMAT1A);
   		DO j=1 to nrow(IMPMAT1A);
      			IMPMAT1A[i,j] = IMPMAT1B[ct];
      			ct = ct + 1;
   		END;
	END;
END;
ELSE DO;
	PRINT "The first model implied covariance matrix entered does not conform to the # of variables specified";
END;

/*Make the implied matrix a correlation matrix; */
tempdiag = J(NROW(IMPMAT1A),NCOL(IMPMAT1A),0);
DO i = 1 to NROW(IMPMAT1A);
	DO j = 1 to NCOL(IMPMAT1A);
		IF i = j THEN tempdiag[i,j]=1/sqrt(IMPMAT1A[i,j]);
	END;
END;
IMPMAT = tempdiag*IMPMAT1A*tempdiag;

/* Now adjust the magnitude of this matrix; */
DO i = 1 to NROW(IMPMAT);
	DO j = 1 to NCOL(IMPMAT);
		IMPMAT[i,j]=IMPMAT[i,j]*.005;
	END;
END;
FREE tempmat tempdiag IMPMAT1A IMPMAT1B;


/* Model implied cov matrix for the model with FEWER vanishing tetrads;
Convert this vector into a square matrix to use in CTA;
Symmetrize this matrix, if not already symmetric; */

use imp2;
read all into IMPMAT2B;

IF NCOL(IMPMAT2B) = (&vars*(&vars+1)/2) THEN DO;
/*IF lowdiag = 1 THEN DO; */
	IMPMAT2A=J(&vars,&vars,0);
	ct = 1;
	DO i=1 to nrow(IMPMAT2A);
   		DO j=1 to nrow(IMPMAT2A);
       		IF i < j then DO;
       			IMPMAT2A[i,j] = 0;
       		END;
   			ELSE DO;
       			IMPMAT2A[i,j] = IMPMAT2B[ct];
       			ct = ct + 1;
       		END;
   		END;
	END;
	tempmat=IMPMAT2A+IMPMAT2A`;
	DO i=1 to nrow(tempmat);
   		DO j=1 to nrow(tempmat);
			IF i=j then tempmat[i,j]= tempmat[i,j]/ 2;
		END;
	END;
	IMPMAT2A = tempmat;
END;
ELSE IF NCOL(IMPMAT2B) = &vars * &vars THEN DO;
	/*Convert this vector into a square matrix to use in CTA; */
	IMPMAT2A=J(&vars,&vars,0);
	ct = 1;
	DO i=1 to nrow(IMPMAT2A);
   		DO j=1 to nrow(IMPMAT2A);
      			IMPMAT2A[i,j] = IMPMAT2B[ct];
      			ct = ct + 1;
   		END;
	END;
END;
ELSE DO;
	PRINT "The second model implied covariance matrix entered does not conform to the # of variables specified";
END;

/*Make the implied matrix a correlation matrix; */
tempdiag = J(NROW(IMPMAT2A),NCOL(IMPMAT2A),0);
DO i = 1 to NROW(IMPMAT2A);
	DO j = 1 to NCOL(IMPMAT2A);
		IF i = j THEN tempdiag[i,j]=1/sqrt(IMPMAT2A[i,j]);
	END;
END;
IMPMAT2 = tempdiag*IMPMAT2A*tempdiag;

/* Now adjust the magnitude of this matrix; */
DO i = 1 to NROW(IMPMAT2);
	DO j = 1 to NCOL(IMPMAT2);
		IMPMAT2[i,j]=IMPMAT2[i,j]*.005;
	END;
END;
FREE IMPMAT2B IMPMAT2a tempmat tempdiag;

CORR = 0;
IF SUM(DIAG(SAMPMAT))=NCOL(SAMPMAT) THEN CORR = 1;

k = NCOL(SAMPMAT);


/* computing sample tetrad residuals */
tetres = {0}; tetlabel = {0 0 0 0 0 0 0 0};
DO a = 1 to k;
	DO b = a+1 to k;
		DO c = b+1 to k;
			DO d = c+1 to k;
			/* NOTE: the 3rd and 4th elements of the 2nd and 3rd tetrads should be reversed to 
			allow calculating the SE's correctly down below--but then would need to sort
			the elements within the tetrads at a later point; */
			tetres = tetres//(SAMPMAT[a,b]*SAMPMAT[c,d] - SAMPMAT[a,c]*SAMPMAT[b,d]);
			tetlabel = tetlabel//(a||b||c||d||a||c||b||d);
			tetres = tetres//(SAMPMAT[a,c]*SAMPMAT[b,d] - SAMPMAT[a,d]*SAMPMAT[b,c]);
			tetlabel = tetlabel//(a||c||b||d||a||d||b||c);
			tetres = tetres//(SAMPMAT[a,b]*SAMPMAT[c,d] - SAMPMAT[a,d]*SAMPMAT[b,c]);
			tetlabel = tetlabel//(a||b||c||d||a||d||b||c);
			END;
		END;
	END;
END;
tetres = tetres[2:NROW(tetres),];
tetlabel = tetlabel[2:NROW(tetlabel),];

/* computing implied tetrad residuals */
tetresI = {0};
DO a = 1 to k;
	DO b = a+1 to k;
		DO c = b+1 to k;
			DO d = c+1 to k;
 			tetresI = tetresI//(IMPMAT[a,b]*IMPMAT[c,d] - IMPMAT[a,c]*IMPMAT[b,d]);
			tetresI = tetresI//(IMPMAT[a,c]*IMPMAT[b,d] - IMPMAT[a,d]*IMPMAT[b,c]);
			tetresI = tetresI//(IMPMAT[a,b]*IMPMAT[c,d] - IMPMAT[a,d]*IMPMAT[b,c]);
			END;
		END;
	END;
END;
tetresI = tetresI[2:NROW(tetresI),];

/* computing implied tetrad residuals for RESTRICTED model*/
tetresI2 = {0};
DO a = 1 to k;
	DO b = a+1 to k;
		DO c = b+1 to k;
			DO d = c+1 to k;
 			tetresI2 = tetresI2//(IMPMAT2[a,b]*IMPMAT2[c,d] - IMPMAT2[a,c]*IMPMAT2[b,d]);
			tetresI2 = tetresI2//(IMPMAT2[a,c]*IMPMAT2[b,d] - IMPMAT2[a,d]*IMPMAT2[b,c]);
			tetresI2 = tetresI2//(IMPMAT2[a,b]*IMPMAT2[c,d] - IMPMAT2[a,d]*IMPMAT2[b,c]);
			END;
		END;
	END;
END;
tetresI2 = tetresI2[2:NROW(tetresI2),];

/* This is used when bringing in an external asymptotic cov matrix; */
acmlabel = J(NROW(tetlabel),.5*(NCOL(tetlabel)),'aaaa');
i = 1;
DO i = 1 to NROW(tetlabel);
	k = 0;
	DO j = 1 to 7 BY 2;
		k = k+1;
		IF tetlabel[i,j] < 10 & tetlabel[i,j+1] < 10 THEN DO;
			acmlabel[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
		END;
		ELSE DO;
			IF tetlabel[i,j] >= 10 & tetlabel[i,j+1]>=10 THEN
				acmlabel[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
			ELSE DO;
				IF tetlabel[i,j] >= 10 THEN
					acmlabel[i,k] = ROWCATC('0'||CHAR(tetlabel[i,j+1])||CHAR(tetlabel[i,j]));
				ELSE acmlabel[i,k] = ROWCATC(CHAR(tetlabel[i,j+1])||'0'||CHAR(tetlabel[i,j]));
			END;
		END;
	END;
END;
acmcov = UNIQUE(SHAPE(acmlabel,1)); *row vector of unique covariances;
FREE acmlabel;

/* find model-implied vanishing tetrads & calculate t-values */
i = 1;
DO i = 1 to NROW(tetlabel);
 a = tetlabel[i,1]; b = tetlabel[i,2]; c = tetlabel[i,3]; d = tetlabel[i,4];
 IF MOD(i,3) = 2 THEN DO;
 	temp = c; c = d; d = temp;
 END;
 IF MOD(i,3) = 0 THEN DO;
 	temp = c; c = d; d = temp;
 END;
 AVari = 1/&N*(Impmat[b,d]**2 * Impmat[a,a] * Impmat[c,c] + 
              Impmat[a,c]**2 * Impmat[d,d] * Impmat[b,b] + 
			  Impmat[c,d]**2 * Impmat[a,a] * Impmat[b,b] +
              Impmat[a,b]**2 * Impmat[d,d] * Impmat[c,c] +
              2*(Impmat[b,d]*Impmat[a,c]*Impmat[a,d]*Impmat[b,c] -
                 Impmat[b,d]*Impmat[c,d]*Impmat[a,a]*Impmat[b,c] -
                 Impmat[b,d]*Impmat[a,b]*Impmat[a,d]*Impmat[c,c] -
                 Impmat[a,c]*Impmat[c,d]*Impmat[a,d]*Impmat[b,b] -
                 Impmat[a,c]*Impmat[a,b]*Impmat[d,d]*Impmat[b,c] +
                 Impmat[c,d]*Impmat[a,b]*Impmat[a,d]*Impmat[b,c]) +
              2*((Impmat[b,d]*Impmat[a,c]-Impmat[a,b]*Impmat[c,d])**2));
 IF i = 1 THEN AVAR = Avari;
 ELSE AVar = Avar//Avari;
END;
FREE Avari;

/* find model-implied vanishing tetrads for model with FEWER vanishing tetrads 
		& calculate t-values */
i = 1;
DO i = 1 to NROW(tetlabel);
 a = tetlabel[i,1]; b = tetlabel[i,2]; c = tetlabel[i,3]; d = tetlabel[i,4];
 IF MOD(i,3) = 2 THEN DO;
 	temp = c; c = d; d = temp;
 END;
 IF MOD(i,3) = 0 THEN DO;
 	temp = c; c = d; d = temp;
 END;
 AVari2 = 1/&N*(Impmat2[b,d]**2 * Impmat2[a,a] * Impmat2[c,c] + 
              Impmat2[a,c]**2 * Impmat2[d,d] * Impmat2[b,b] + 
			  Impmat2[c,d]**2 * Impmat2[a,a] * Impmat2[b,b] +
              Impmat2[a,b]**2 * Impmat2[d,d] * Impmat2[c,c] +
              2*(Impmat2[b,d]*Impmat2[a,c]*Impmat2[a,d]*Impmat2[b,c] -
                 Impmat2[b,d]*Impmat2[c,d]*Impmat2[a,a]*Impmat2[b,c] -
                 Impmat2[b,d]*Impmat2[a,b]*Impmat2[a,d]*Impmat2[c,c] -
                 Impmat2[a,c]*Impmat2[c,d]*Impmat2[a,d]*Impmat2[b,b] -
                 Impmat2[a,c]*Impmat2[a,b]*Impmat2[d,d]*Impmat2[b,c] +
                 Impmat2[c,d]*Impmat2[a,b]*Impmat2[a,d]*Impmat2[b,c]) +
              2*((Impmat2[b,d]*Impmat2[a,c]-Impmat2[a,b]*Impmat2[c,d])**2));
 IF i = 1 THEN AVAR2 = Avari2;
 ELSE AVar2 = Avar2//Avari2;
END;
FREE AVari2;


/* This section is used for nested tetrad tests; */
IF &nesttest = 1 THEN DO;
	/* Check whether the models are tetrad nested, and which is bigger; */
	temp1 = J(1,1,0);
	temp2 = J(1,1,0);
	temp3 = J(1,1,0);
	difftet = J(1,NCOL(tetlabel)+2,0);
	badtet = J(1,NCOL(tetlabel)+2,0);

	tvalueI = J(NROW(tetresI),3,0);
	DO i = 1 to NROW(tvalueI);
		tvalueI[i,1]= tetresI[i,1];
		tvalueI[i,2]= AVar[i,1];
		tvalueI[i,3]=tvalueI[i,1]/sqrt(tvalueI[i,2]);
	END;

	tvalueI2 = J(NROW(tetresI2),3,0);
	DO i = 1 to NROW(tvalueI2);
		tvalueI2[i,1]= tetresI2[i,1];
		tvalueI2[i,2]= AVar2[i,1];
		tvalueI2[i,3]=tvalueI2[i,1]/sqrt(tvalueI2[i,2]);
	END;
	DO i = 1 to NROW(tetlabel);
		/* This is the t-value criterion for defining as vanishing (keep non-vanishing);
		This loop finds vanishing tetrads in the first model but not in the second one; */

		IF ABS(tvalueI[i,3]) < 0.001 & ABS(tvalueI2[i,3]) >= 0.001 THEN DO;
			temp1 = temp1 || i;
			temptet = tetlabel[i,]||tvalueI[i,3]||tvalueI2[i,3];
			difftet = difftet // temptet;
		END;
		/* This loop finds vanishing tetrads in the second model but not in the first one; */
		IF ABS(tvalueI[i,3]) >= 0.001 & ABS(tvalueI2[i,3]) < 0.001 THEN DO;
			temp2 = temp2 || i;
			temptet2 = tetlabel[i,]||tvalueI[i,3]||tvalueI2[i,3];
			badtet = badtet // temptet2;
		END;
		/* This loop finds vanishing tetrads implied by both models; */
		IF ABS(tvalueI[i,3]) < 0.001 & ABS(tvalueI2[i,3]) < 0.001 THEN DO;
			temp3 = temp3 || i;
		END;
	END;

	/* **Find ERROR in specification (not nested); */
	notnest = 0;
	if NCOL(temp1) >= 2 & NCOL(temp2) >=2 then do;
		notnest = 1;
	END;

	if NCOL(temp1) = 1 & NCOL(temp2) = 1 & NCOL(temp3) >=2 then do;
		notnest = 2;
	END;

	if NCOL(temp1) = 1 & NCOL(temp3) = 1 then do;
		notnest = 3;
		zerotet = 1;
	END;

	if NCOL(temp2) = 1 & NCOL(temp3) = 1 then do;
		notnest = 4;
		zerotet = 1;
	END;

	/* These are models in which the order of the model implied cov's was reversed between smaller 
	and larger models, so need to flip them here to correct operator error; */

	matswitch = 0;
	if NCOL(temp1) = 1 & NCOL(temp2) >=2 then do;
		matswitch = 1;
		temp21 = AVar2;
		temp22 = tetresI2;
		temp23 = IMPMAT2;
		AVar2 = AVar;
		tetresI2 = tetresI;
		IMPMAT2 = IMPMAT;
		AVar = temp21;
		tetresI = temp22;
		IMPMAT = temp23;
	END;
END;
FREE temp1 temp2 temp3 difftet badtet temp21 temp22 temp23;



/* Combining the residuals with the implied tetrads and computing t-values; */
tvalueI = J(NROW(tetresI),6,0);
DO i = 1 to NROW(tvalueI);
	tvalueI[i,1]= tetresI[i,1];
	tvalueI[i,2]= AVar[i,1];
	tvalueI[i,3]=tvalueI[i,1]/sqrt(tvalueI[i,2]);
	tvalueI[i,4]= tetresI2[i,1];
	tvalueI[i,5]= AVar2[i,1];
	tvalueI[i,6]=tvalueI[i,4]/sqrt(tvalueI[i,5]);
END;

/* Computing the variances for the sample tetrads; */
DO i = 1 to NROW(tetlabel);
 a = tetlabel[i,1]; b = tetlabel[i,2]; c = tetlabel[i,3]; d = tetlabel[i,4];
 IF MOD(i,3) = 2 THEN DO;
 	temp = c; c = d; d = temp;
 END;
 IF MOD(i,3) = 0 THEN DO;
 	temp = c; c = d; d = temp;
 END;

         AVars = (1/N) *
               (Sampmat[d,b]*Sampmat[d,b]*Sampmat[a,a]*Sampmat[c,c] +
                Sampmat[a,c]*Sampmat[a,c]*Sampmat[d,d]*Sampmat[b,b] +
                Sampmat[d,c]*Sampmat[d,c]*Sampmat[a,a]*Sampmat[b,b] +
                Sampmat[a,b]*Sampmat[a,b]*Sampmat[d,d]*Sampmat[c,c] +
            2 *(Sampmat[d,b]*Sampmat[a,c]*Sampmat[a,d]*Sampmat[b,c] -
                Sampmat[d,b]*Sampmat[d,c]*Sampmat[a,a]*Sampmat[b,c] -
                Sampmat[d,b]*Sampmat[a,b]*Sampmat[a,d]*Sampmat[c,c] -
                Sampmat[a,c]*Sampmat[d,c]*Sampmat[d,a]*Sampmat[b,b] -
                Sampmat[a,c]*Sampmat[a,b]*Sampmat[d,d]*Sampmat[b,c] +
                Sampmat[d,c]*Sampmat[a,b]*Sampmat[a,d]*Sampmat[b,c])+
           2 *((Sampmat[d,b]*Sampmat[a,c]-Sampmat[a,b]*Sampmat[d,c])**2));

IF i = 1 THEN AVart = Avars;
 ELSE AVart = Avart//Avars;
END;
FREE AVars;
/* Combining the residuals with the tetrads and computing t-values; */
tvalue = J(NROW(tetres),3,0);
do i = 1 to NROW(tvalue);
	tvalue[i,1]= tetres[i,1];
	tvalue[i,2]= AVart[i,1];
	tvalue[i,3]=tvalue[i,1]/sqrt(tvalue[i,2]);
END;

/* Now keep the tetrads that are implied vanishing by the model; */
tetlabel = tetlabel||tvalueI;
tetres = tetres||tvalueI;
temp1 = J(1,8,0);
temp2 = J(1,1,0);
temp3 = J(1,8,0);
temp4 = J(1,1,0);
temp5 = J(1,8,0);
temp6 = J(1,1,0);

DO i = 1 to NROW(tetlabel);
	/* This is the t-value criterion for defining as vanishing (keep non-vanishing);
	This loop finds all vanishing tetrads; */
	IF ABS(tetlabel[i,11]) < 0.001 THEN DO;
		temp1 = temp1 // tetlabel[i,1:8];
		temp2 = temp2 // tetres[i,1];
	END;
	/* This loop finds vanishing tetrads in the larger model but not in the smaller one; */
	IF ABS(tetlabel[i,11]) < 0.001 & ABS(tetlabel[i,14]) >= 0.001 THEN DO;
		temp3 = temp3 // tetlabel[i,1:8];
		temp4 = temp4 // tetres[i,1];
	END;
	/* This loop finds vanishing tetrads implied by both models; */
	IF ABS(tetlabel[i,11]) < 0.001 & ABS(tetlabel[i,14]) < 0.001 THEN DO;
		temp5 = temp5 // tetlabel[i,1:8];
		temp6 = temp6 // tetres[i,1];
	END;
END;

/*Comparing nested models--- Check that they are nested; */
notnest2 = 0;
IF &nesttest = 1 & NROW(temp3) = 1 THEN DO;
	notnest2 = 1;
END;
IF NROW(temp1) = 1 THEN DO;
	print "There are no vanishing tetrads in either model";
	/* FIX THIS */
	/* NEED TO JUMP OUT OF THE PROGRAM AT THIS POINT */
END;
ELSE IF &nesttest = 1 THEN DO;
	temp3 = temp3[2:NROW(temp3),];
	temp4 = temp4[2:NROW(temp4),];
END;
ELSE DO;
	/* This is the non-nested test: so keep all of its tetrads (temp1 and temp2); */
	temp3 = temp1[2:NROW(temp1),];
	temp4 = temp2[2:NROW(temp2),];
END;
IF &nesttest = 1 & zerotet = 0 THEN DO;
	tetlabel = temp5[2:NROW(temp5),];
	tetres = temp6[2:NROW(temp6),];
END;
/* This is a nested test where one model has NO vanishing tetrads; */
ELSE IF &nesttest = 1 & zerotet = 1 THEN DO;
	tetlabel = temp3;
	tetres = temp4;
END;
ELSE DO;
	/* This is a non-nested test; */
	tetlabel_nonnest = temp3;
	tetres_nonnest = temp4;
	initlabel_nonnest = tetlabel_nonnest;
	initres_nonnest = tetres_nonnest; 
END;
/* Don't free up temp3 temp4 yet, as they are needed below; */
FREE temp1 temp2 temp5 temp6;
	if notnest2 = 1 then print "These models are NOT nested---same # of df";
	if matswitch = 0 then print "Model implied matrices are in correct order";
	if matswitch = 1 then print "NOTE: Model implied matrices were incorrectly entered for large and small models---they have been reversed";
	if notnest = 0 then print "Models are tetrad nested";
	if notnest = 1 then print "ERROR:  These models are NOT tetrad nested.  Nested test is INAPPROPRIATE";
	if notnest = 2 then print "ERROR:  These models contain the same tetrads.  Nested test is INAPPROPRIATE";
	if notnest = 3 then print "NOTE:  The first model contains no tetrads; tetrad chi-square for first model also represents nested test";
	if notnest = 4 then print "NOTE:  The second model contains no tetrads; tetrad chi-square for first model also represents nested test";

initlabel = tetlabel;
initres = tetres;

DO r = 1 to &reps;
	RUN TETRAD;
END;


/* Fixing the formatting of the output */
IF &reps > 1 THEN DO;
	resultsch = J(NROW(results),NCOL(results),'aaaaaaaaa');
	IF &nesttest = 1 THEN DO;
		DO q = 1 to NROW(results);
			resultsch[q,] = CHAR(results[q,1],4,0)||CHAR(results[q,2],8,2)||CHAR(results[q,3],4,0)||
			CHAR(results[q,4],8,4)||CHAR(results[q,5],8,2)||CHAR(results[q,6],4,0)||CHAR(results[q,7],8,4)
			||CHAR(results[q,8],8,2)||CHAR(results[q,9],4,0)||CHAR(results[q,10],8,4);
		END;
  	END;
  	ELSE DO;
		DO q = 1 to NROW(results);
			resultsch[q,] = CHAR(results[q,1],4,0)||CHAR(results[q,2],8,2)||CHAR(results[q,3],4,0)||
			CHAR(results[q,4],8,4);
		END;
	END;
END;


IF &reps > 1 THEN DO;
	PRINT Resultsch[label='results' colname=clabel];
	CREATE Outfile FROM RESULTS[colname=clabel];
		APPEND FROM RESULTS;
END;

FINISH;
RUN;
QUIT;

%mend ctanest1;