How to use JSL to save T² from PCA platform (or how to calculate it)

SDF1

Hi All,

 

  I am trying to get some JSL code working to save T² from the Principal Component platform, like how you can do it by calling <<Save DModX(n), where n is the number of PCs to use.

 

  I know that the call << Save T Square will work with multivariate control chart, but it doesn't seem to work with the PCA or Multivariate platforms. I need to use the PCA platform to generate the T² formula because in the PCA platform, you can actually specify the number of PCs to use in the calculation, it's not the normal T² statistic, like the Mahalanobis distance squared. Even the call << T²(1, Save T Square) sent to the report window doesn't work.

 

1. Is it possible to use scripts to specify the number of PCs to use and save the T² formula to the data table? So far, I have not found an option for this.

 

2. If not, how to best calculate it? I know the calculation is defined here, but the formula that's saved to data table using the red hot-button seems a bit different.

 

  According to the saved formula (manual save) T² = Vec Quadratic(S, X)= Diag(X*S*X`), where X is the standardized vector of data, X` is its transpose, and S is the inverse covariance matrix. This is the standard T², the square of the Mahalanobis distance.

 

  However, according to the documentation for the ith T² value from the PCA platform, T_i² = [X_ci*P_a]*L-1*[P`_a*X`_ci], where X_ci is the ith standardized observation, P_a is a matrix containing the first "a" eigenvalues (i.e. possibly not all of them), and L-1 is a diagonal matrix containing the first A eigenvalues. In principle, this should be equal to Vec Quadratic(L-1, [X_ci*P_a]), except instead of using the inverse covariance matrix, S, it's using the inverse axa eigenvalue matrix of the first "a" PCs, and X is now X_ci*P_a]. Or so, it would appear.

 

  It seems rather complicated and unnecessary to re-engineer this problem in order to save the n-principal component T² when JMP already calculates all this for you. I just want to automate the saving process, but I can't seem to get option 1 to work, which means option 2 is the next best bet. Even then, trying to recover the formula in order to make a New Column with a formula in it seems a bit much.

 

  I'm attaching a modified data table of the Iris.jmp file that shows the saved T² formulas from the multivariate platform (and how T² = Mahal.Dist.²), and then T² from the PCA platform where the formula is determined using all 4 PCs and then only to PCs. As you can see by looking at the formulas, the S matrix in the Vec Quadratic(S,X) is different for the T²s calculated using 2 PCs vs 4 PCs.

 

  If anyone has insight into how to script getting the T² with a different number of PCs or how to properly formulate the column formula to calculate the T², I would be grateful for that insight.

 

Thanks for your feedback and thoughts!,

DS

 

Iris_mod.jmp

SDF1

Hi @Mark_Bailey ,

 

  Thanks for your feedback and including the JSL from the developer of the platform. Well, as it turns out, the solution to what I was after (option 1) is actually possible with JSL!! And so much simpler, too. I just needed to see it in the code that you sent.  The critical part is here:

 

Outlier Analysis( 1, Number of Components( N ), Save T Square ),

  By changing "N" in the Number of Components input, and adding the Save T Square, it will actually save the T² formula to the data table for the N components that you define. I didn't know that within the Outlier Analysis() call, you could send it the request to "Save T Square". Doing that along with the N components does exactly what I was hoping to do, thanks for the help! No need to put in a Wish List or anything, JMP can already do it for you!

 

 

  I'm attaching the Iris_mod.jmp data table with scripts to that calculate the T² value for 2, 3, and 4 PCs. It is good to know how to efficiently calculate the T² value should I ever need to. The code I was just writing now to do the calculation was rather long and cumbersome, although I believe it would have worked once I got everything straightened out. However, I prefer the above solution -- it's only a few extra characters of code!

 

  Upon closer inspection of the code you sent, I have noticed that the T² calculated from the PCA platform and the "manual" calculation in the code are different for the same number of PCs. Any clue as to why? Either the "manual" calculation is "wrong" or performed differently than how the PCA does it, or something else is going on. I'm also attaching another data table Iris_mod_different_T2.jmp that shows the T² results for the same number of PCs is different from the PCA platform and the "manual" calculation. Until further notice, I prefer to use the T² formula from the PCA platform. I would like to know if there something different about the "manual" calculation and why the values are not the same.

 

Thanks again!,

DS

 

 

 

P.S. If anyone is interested, the code I've been working on is below. It is an interactive code that allows the user to select the columns to perform the PCA on and then dynamically select the number of PCs to use for extracting out the DModX and T² formulas. You can even save the selected number of PCs to the data table if you want for other analyses.

 

  A big thanks to @jthi for help with some other aspects of this code in order to get things to work as intended, like <<Lower() for the Plot Col Box() and <<Set Selectable Rows(1) for the Table Box() display boxes. The only thing I couldn't figure out is how to limit the selection of rows in the Table Box() to just one. In the Col List Box() environment, you can set the Max Items(1) to 1 and this only allows the person to select one item at a time. Any hints on how to do that for the code below is much appreciated. Right now, it just takes the max PC number if the user selects more than one row in the Table Box().

Names Default To Here( 1 );

//Clear Symbols();

//Function to choose the data table to model
choose_data_table = Function( {},
	list = {};
	For( i = 1, i <= N Table(), i++,
		list[i] = Data Table( i ) << get name
	);
	win = New Window( "Select a data table",
		<<Modal,
		hb = H List Box(
			Panel Box( "Choose a data table to extract T² and DModX",
				dt = List Box( list, max selected( 1 ), dtchosen = Data Table( (dt << get selected)[1] ) )
			)

		)
	);
);

// Expression to clear all current settings
clearRoles = Expr( ColListY << RemoveAll );

// Expression to store the current settings in global variables
recallRoles = Expr(
	::ycolRecall = ColListY << GetItems;
	::dtchosenRecall = dtchosen;
);

//Expression to get T² and DModX after selecting the corresponding PC number
Get_T_DModX = Expr(

	T2_text = "T² for " || Char( N_PCs ) || " Principal Components";

	PCRptOutlier = dtchosen << Principal Components(
		Y( Eval( tpListY ) ),
		Estimation Method( "Default" ),
		"on Correlations",
		Eigenvectors( 1 ),
		Eigenvalues( 1 ),
		Outlier Analysis( 1, Number of Components( N_PCs ), Save T Square ),
		SendToReport(
			Dispatch( {}, "Summary Plots", OutlineBox, {Close( 1 )} ),
			Dispatch( {}, "Eigenvalues", OutlineBox, {Close( 1 )} ),
			Dispatch( {"Outlier Analysis"}, "PrinComp", FrameBox, {Reference Line Order( 2 ), Grid Line Order( 1 )} )
		),
		Invisible
	);<BR />    PCRptOutlier << Save DModX( N_PCs );
	If( savePCs == 1,
		PCRptOutlier << Save Principal Components( N_PCs )
	);	
	PCRptOutlier<<CloseWindow;
);

//Expressoin to get PCs
Metrics = Expr(
	PCRpt = dtchosen << Principal Components(
		Y( Eval( tpListY ) ),
		Estimation Method( "Default" ),
		"on Correlations",
		Eigenvalues( 1 ),
		Invisible
	);
	
	PC_dt = Report( PCRpt )[Outline Box( "Principal Components: on Correlations" )][Outline Box( "Eigenvalues" )][
	Table Box( 1 )] << Make Into Data Table;
	
	PC_dt << Set Name( "PC Table Results" );
	PC_n = [];
	PC_eigv = [];
	PC_per = [];
	PC_cump = [];
	For( i = 1, i <= N Rows( PC_dt ), i++,
		Insert Into( PC_n, PC_dt:"Number"n[i] );
		Insert Into( PC_eigv, Round( PC_dt:"Eigenvalue"n[i], 5 ) );
		Insert Into( PC_per, Round( PC_dt:"Percent"n[i], 5 ) );
		Insert Into( PC_cump, Round( PC_dt:"cumulative Percent"n[i], 5 ) );
	);
	
	PC_Rpt_win = New Window( "PC Results",
		<<Return Result,
		<<On Validate,
		Outline Box( "Eigenvalues & cumulative Percent",
			H List Box(
				V List Box(
					tb = Table Box(
						Number Col Box( "PC Number", Eval List( PC_n ) ),
						Number Col Box( "Eigenvalue", Eval List( PC_eigv ), <<Set Format( "Fixed Dec", 6, 6 ) ),
						Number Col Box( "Percent", Eval List( PC_per ), <<Set Format( "Fixed Dec", 6, 6 ) ),
						Plot Col Box( "",
							Eval List( PC_per ),
							<<Set Scale( 0, 100, "Percent" ),
							<<Lower( PC_cump ),
							<<Upper( PC_cump )
						),
						Number Col Box( "cumulative Percent", Eval List( PC_cump ), <<Set Format( "Fixed Dec", 6, 6 ) ),
						<<Set Shade Headings( 1 ),
						<<Set Heading Column Borders( 1 ),
						<<Set Selectable Rows( 1 ), 

					),
					H List Box( Spacer Box( Size( 58, 0 ) ), savePCs_input = Check Box( "Save PCs?" ) )
				),
				Panel Box( "Action",
					Lineup Box( N Col( 1 ),
						Text Box( "Get T² and DModX" ),
						Button Box( "OK",
							savePCs = PC_Rpt_win[CheckBoxBox( 1 )] << Get;
							Get_T_DModX;
							PC_Rpt_win << Close Window;
						),
						Button Box( "Cancel", PC_Rpt_win << Close Window ),
						Spacer Box( Size( 0, 10 ) ),
						Button Box( "Relaunch",
							PC_Rpt_win << Close Window;
							Firstwin;
						),
						Spacer Box( Size( 0, 10 ) ),
						Button Box( "Help", Web( "https://www.jmp.com/en_ch/support/online-help-search.html?q=*%3A*" ) )
					)
				)
			)
		)
	);
	Close( PC_dt, No Save );
	PCRpt << CloseWindow;
	tb << Set Row change Function( Function( {this}, N_PCs = Maximum( this << Get Selected Rows ) ) );
);

//Dialog Window.
ModMetrics = Expr(
	DistMetricswin = New Window( "Extracting Model Eval Metrics, T² and DModX",
		<<Return Result,
		<<On Validate,
		Border Box( Left( 3 ), top( 2 ),
			Outline Box( "GUI to simplify extracting T² and DModX",
				<<Set Font Size( 12 ),
				H List Box(
					V List Box( Panel Box( "Select Factors", colListData = Col List Box( all, Grouped, nLines( 12 ) ) ) ),
					Panel Box( "Cast columns into their role",
						V List Box(
							Lineup Box( N Col( 2 ), Spacing( 3 ),
								Button Box( "Variables", colListY << Append( colListData << Get Selected ) ),
								colListY = Col List Box( nLines( 11 ), "numeric", Min Items( 2 ) )
							),
							Spacer Box( Size( 0, 11 ) )
						)
					),
					Panel Box( "Action",
						Lineup Box( N Col( 1 ),
							Button Box( "OK",
								recallRoles;
								tpListY = colListY << Get Items;
								Metrics;
								DistMetricswin << Close Window;
							),
							Button Box( "Cancel", DistMetricswin << Close Window ),
							Spacer Box( Size( 0, 26 ) ),
							Button Box( "Remove", tpListY << Remove Selected ),
							Button Box( "Recall",
								clearRoles;
								Try( colListY << Append( ::ycolRecall ) );
							),
							Button Box( "Relaunch",
								DistMetricswin << Close Window;
								Firstwin;
							),
							Spacer Box( Size( 0, 26 ) ),
							Button Box( "Help", Web( "https://www.jmp.com/en_ch/support/online-help-search.html?q=*%3A*" ) )
						)
					), 

				)
			)
		)
	)
);


//Interactive dialogue window to start Saving T² and DModX
Firstwin = Expr(
	dtlist = Get Data Table List();
	If(
		N Items( dtlist ) == 0,
			Try( dt_O = Open(), Throw( "No data found" ) );
			dtchosen = dt_O;,
		N Items( dtlist ) > 0, Try( choose_data_table, Throw( "No data found" ) )
	);
	ModMetrics;
);

Firstwin;

Iris_mod.jmp

Iris_mod_different_T2.jmp