In this chapter, we delve into the sophisticated realm of multivariate statistical methods, focusing on Principal Component Analysis (PCA) and Projection to Latent Structures (PLS), as pivotal tools for unraveling the complexity of process data analytics . By anchoring these statistical techniques within the framework of polyethylene manufacturing processes, we aim to illuminate their exceptional utility and novelty in addressing the multifaceted challenges inherent in process optimiz ation and quality control. The discourse begins by introducing PCA, not merely as a statistical tool, but as a fundamental cornerstone for the analytical examination of process variables. Through a meticulously designed workshop, we demonstrate the applica tion of PCA in dissecting the intricate web of variables influencing the quality and conversion rates of Low -Density Polyethylene (LDPE) production in a two -zone tubular reactor. The integration of Aspen ProMV as a practical tool for PCA applications exemp lifies the seamless bridge between statistical theory and industrial application, emphasizing the method's accessibility and relevance to both academia and industry. Transitioning to PLS, the chapter articulates its differentiation from PCA by its ability to simultaneously handle datasets comprising both process variables (X) and product quality variables (Y), offering a holistic view of the manufacturing process. Through pragmatic workshops, we showcase PLS's robustness in application to challenges such as melt index prediction and causal analysis in High -Density Polyethylene (HDPE) manufacturing, underscoring its adaptability to complex industrial datasets, including those with measurement time lags. The exploration extends to the nuanced application of th ese multivariate statistical methods to batch polymer processes. Here, we introduce a novel batch -wise unfolding approach via multiway PCA and PLS, expanding the frontier of statistical applications in process data analytics. This chapter transcends the conventional boundaries of statistical applications, highlighting the transformative impact of PCA and PLS in the domain of process data analytics. It aspires to foster a deeper understanding a nd appreciation of these statistical methods, encouraging their broader adoption and adaptation in optimizing manufacturing processes and enhancing product quality. This contribution not only reaffirms the critical role of advanced statistical techniques i n the scientific community but also underscores their practical significance in improving industrial operations and outcomes. This is a preprint version of a chapter from our book - Integrated Process Modeling, Advanced Control and Data Analytics for Optimizing Polyolefin Manufacturing . Please cite the original work if referenced [31,32 ]. 9.1 Introduction to Principal Component Analysis (PCA) This Chapter 9 focuses on the use of multivariate statistics . Sections 9.1 introduces an important multivariate statistics tool in process data analytics, namely, principal component analysis (PCA). Section 9.2 presen ts a hands -on workshop on the application of PCA for analyzing the process variables that affect quality and conversion of LDPE product from a two -zone tubular reactor. We introduce the use of the software tool, Aspen ProMV, for multivariate statistics app lications, available to universities at low cost. Section 9.3 introduces the projection to latent structures or partial least squares (PLS) . A key difference between PCA and PLS is that PCA involves only datasets of process variables (X) or deals with the X-space; PLS involves datasets of both process variables (X) and product quality variables (Y), or deals with both X -space and Y -space. Section 9.4 presents two hands -on workshops of applying PLS to the LDPE problem of Section 9.2, and to the melt index pr ediction and causal analysis of a HDPE manufacturing process. Section 9.5 introduces PLS for process data analytics with measurement time lags and includes a hands -on workshop of PLS for a HDPE process for the melt index prediction and causal analysis, inc luding the effect of time lags on melt index measurements. Section 9.6 covers the process data analytics for batch polymer processes and presents a hands -on workshop to demonstrate the multiway PCA and PLS methodology, particularly the batch -wise unfolding approach, for data analytics. Section 9.7 gives the bibliography and suggested further reading Beginning in late 1980 to early 1990, chemical engineers have been paying an increasing attention to th e emerging topics of artificial intelligence, neural computing, multivariate statistics, machine learning and big data analytics, and their applications to bioprocessing and chemical engineering [1 to 5]. MacGregor and others have demonstrated the signific ant applications of multivariate statistics and big data analytics to optimizing the manufacturing of LDPE, HDPE, Nylon 6 and other polymers [6 to 10]. Multivariate statistical analysis [11 to 13] and its implementation using languages such as Python, R or software such as Aspen ProMV, SAS, JMP, etc. find a growing number of applications to polymer manufacturing. These include: (1) data quality deviation analysis; (2) unit yield analysis; (3) production capacity degradation analysis; (4) offline production optimization (discovery and optimization of key variables); (5) online process monitoring and troubleshooting; and (6) batch process variable analysis. This section introduces the principal component analysis (PCA), following the multivariate statistical analysis textbooks of Johnson and Wichern [11], and Rencher and Christensen [12], and the excellent online book of Dunn [13] , which is continually updated. The online reader is allowed “to freely download, share, adapt, commercialize and attribute” some of the book materials, as long as the reader acknowledges that “Portions of this work are the copyright of Kevin Dunn”. That is exactly what we wish to acknowledge here, as we shall use some of the explanations and figures from reference [13] below. Both text books [11,12] include a chapter of matrix algebra relevant to multivariate statistical analysis. Therefore, we have included an Appendix A, Matrix Algebra in Multivariate Data Analysis and Mode l- Predictive Control Using MATLAB and Python at the end of thi s book. This appendix also includes the basic implementation of the relevant matrix operations and principal component analysis in both MATLAB and in Python. 9.1.1 Introduction to Principal Components We follow [14] to illustrate the concept of principal components. Figure 9.1 shows a 3D image of some process data. When projecting the same data onto a 2D plane in Figure 9.2, we are unable to observe the same 3D relationship. However, we can observe sufficient characteristics of the original 3D image in two dimensions if we can identify two linear combinations of process variables x, y and z in order to capture most of the variations in these three process variables. See Figure 9.3. Figure 9.1 The original 3D image of process data . Used with perssion from Aspen Technology, Inc. Figure 9.2 Losing the characteristics of the original 3D image when projecting onto the x -y plane . Used with perssion from Aspen Technology, Inc. Figure 9.3 Retaining the characteristics of the original 3D image when displaying on the two -dimensional plane of latent variables LV1 and LV2 (or principal components 1 and 2) . Used with perssion from Aspen Technology, Inc. In Figure 9.3, we see the characteristics of the original 3D image on a two -dimensional plane of two linear comb inations (LV1 and LV2) of the process variables (x, y and z): LV1 = 0.1658 x + 0.6120 y + 0.7733 z (9.1a) LV2 = -0.9652 x + 0.2615 y (9.1b) We call these linear combinations the latent variables or principal components of the process variables. Principal component analysis (PCA) is a data transformation method that rotates data such that the principal axis of the data is in the direction of maximum variation. See Figure 9.4. We follow the interpretation of [15] here. The first latent variable or first principal component of the process data or observations give n by Eq. (9.1 a) represents the linear combination of the original process variables whose sample variance (see Appendix A, Section A.1.6) is greatest among all possible linear combinations. The second latent variable or second principal component represent s the linear combination of the original process variables that accounts for a maximum proportion of the remaining variance, subject to being uncorrelated with the first principal component. We can define subsequent principal components similarly. Figur e 9.4 An illustration of the principal component that shows the direction of maximum variation of the process data We can view the rotated data on the new principal axes (components). We call the coordinates of the data in this new coordinate system as principal component scores . They are essentially the projection of the data onto the principal axes. As see n in Figure 9.4, the principal components are essentially vectors in the original variable space, and these vectors are called principal component loadings . We quantify both the principal component scores and loadings, and their relationship to the original process data matrix in the following section. 9.1.2 Data Prepro cessing: Mean -Centered and Scaled Process Data Matrix X, Principal Component Score Matrix T, and Principal Component Loading Matrix P We follow [13] for the development of the PCA model. Let us consider a J x K process data matrix X, with K columns of process variables x k (k=1,2 ,…, K), and with each variable x k having J observations or measured values, x 1k, x2k, x3k, …x Jk (or x jk, j = 1, 2,…, J). In Appendix A, Section A.1.7 , we introduce the standardized data matrix, or mean -centered and scaled data matrix Xs, and the correlation coefficient matrix R from the process data matrix X. To correctly carry out PCA, we first preprocess the data. Specifically, we start with a data standardization step to convert the process data matrix X to a standardized data matrix Xs that is mean - centered and scal ed by standard deviation [11,13]. For convenience in eliminating the letter “ s” from a mean -centered and scaled data matrix Xs, we assume in the following discussion that our process data matrix X has already gone through a standardization procedure descri bed in Appendix A, Sections A.1.5 to A.1.7. As we demonstrate in Appendix A, it only takes a single command using Matlab [zscore(X)] or Python [stats.zscore()], to standardize a process data matrix X. 9.1.3 Development of PCA Model We write the standardiz ed data matrix X as a matrix of K process variable vectors: X = [x1 x2 ……..x K] = [𝑥11 ⋯ 𝑥1𝑘 ⋮ ⋱ ⋮ 𝑥𝑗1⋯ 𝑥𝑗𝑘] (9.2) In this matrix, the k -th process variable vector xK is a (J x 1) column vector, [x 1k, x2k, x3k, …x Jk ]’, where J is the numbe r of samples or measurements. The transpose of xK, or xk’, is a (1 x J) observation vector. Figure 9.5 illustrates the projection of the vector xk onto the first principal component vector p1. The score value t k,1 for this observation vector is the distanc e from the origin along the principal component loading vector, p1, to the point where we find the perpendicular projection onto p1 [13]. Figure 9.5 The projection of the kth process variable vect or xk onto the first principal component loading vector p1. tk,1 is the score value of xk on p1. We can write from geometry that: (1) the cosine of an angle in a right -angled triangle is the ratio of the adjacent side to the h ypotenuse; (2) the cosine of the defines the dot product of two vectors . See Eqs. (9.3) and (9.4) : cos θ = (adjacent length)/(hypotenuse) = t k,1 / ‖ xk ‖ (9.3) cos θ = xk’ p1 / ‖ xk ‖ ‖ p1 ‖ (9.4) where ‖ · ‖ represents the length of the enclosed vector, and the length o f the principal component loading vector, ‖ p1 ‖ is 1.0. Therefore, we find: tk,1 = xk’ p1 = x k,1 p1,1 + x k,2 p2,1 + …+ x k,j pj,1 …. + x k,J pJ,1 (9.5) Likewise, we write tk,2 = xk’ p2 = x k,1 p1,2 + x k,2 p2,2 + …+ x k,j pj,2 …. + x k,J pJ,2 (9.6) Generalizing Eq. (9.5) and (9.6), we write the principal component score vector tk resulting from projecting the process data vector xk onto A principal component loading vectors, expressed by the (K x A) loading matr ix P: tk’ = x k’ P (9.7) (1 x A) = (1 x K) (K x A) Lastly, we can represent the projection of the entire process data matrix X in terms of a principal component score matrix T and a principal component loading matrix P: T = X P (9.8) (J x A) = (J x K)(K x A) where J is the number of samples or measurements, A is the number of principal components, and K is the number of process variables. 9.1.4 Prediction Errors from PCA Model Figure 9.6 illustrates the projection of the original data vector xk onto the first principal component vector p1. The best estimate of xk is a vector 𝒙̂𝒌,𝟏along the first principal component loading vector p1 where the original vector is projected. We call this estimate of the data vector, 𝒙̂𝒌,𝟏. We note the distance along the first principal component loading vector p1 is the principal comp onent score t k,1. Based on vector geometry, we represent the error between xk and 𝒙̂𝒌,𝟏as an error vector ek,1. Figure 9.6 The projection of the kth process variable vector xk onto the first principal component vector p1, indicating an es timate of the data vector 𝒙̂𝒌,𝟏, together with an error vector ek,1. tk,1 is the score value of xk on p1. We write the prediction vector as: 𝒙̂𝒌,𝟏′= tk,1 p’1 (9.9) (1 X K) = (1 x 1)(1 x K) and the corresponding prediction error vector is: e’k,1 = x’k - 𝒙̂𝒌,𝟏′ (9.10) (1 x K) = (1 x K) – (1 x K) Adding the second principal component vector p2 , we generalize the prediction vector from Eq. (9.9) as: 𝒙̂𝒌,𝟐′= tk,1 p’1 + tk,2 p’2 (9.11) (1 x K) = (1 x K) + (1 x K) where t k,1 and tk,2 are the score values of xk,2 on p1 and p2, respectively. Extending Eq. (9.11) to A principal component vectors, we write the projector vector of the original data vector xk onto the A principal component loading vectors [ p1 p2 ……p A] or principal component loading matrix P, with tk being the score vector: 𝒙̂𝒌,𝑨′= [tk,1 tk,2…, tk,A ]P’ = t’ k P’ (9.12) (1 x A) = (1 x A) (A x K) We generalize Eq. (9.12) to represent the entire data prediction matrix 𝑿̂ in terms of t he score matrix T and the principal component loading matrix P: 𝑿̂= T P’ (9.13) (J x K ) = (J x A) (A x K) We define the residual vector ek,A for the k -th process variable using A principal components as the difference between the actual and pred icted observations: e’k, A = x’k - 𝒙̂𝒌,𝑨′ = x’k - t’k P’ (9.14) (1 x A) = (1 x A) – (1 x A) Referring to Figure 9.7, we define the row residual or the squared prediction error (SPE) for k -th process variable as: SPE k = ( e’ k, A · ek, A)1/2 = [ (x k,1 – x̂k,1′)2 + (x k,2 – x̂k,2′)2+………+ (x k,A – x̂k,A′)2] 1/2 (9.15) The corresponding vector representation of all SPE k (k= 1 ,2 …K) for all K process variables is SPE = [ SPE 1 SPE 2…, SPE k ]’ (9.16) We write Eq. (9.14) as a prediction error or residual matrix E for all K process variables, J observations per variable, and A principal component loading vectors [ p1 p2 ……p A] or principal component loading matrix P as follows: E = X – 𝑿̂ = X – T P’ or X = T P’ + E (9.17) Figure 9.7 illustrates the relationship between E, X, 𝑿̂, and SPE. Figure 9.7 An illustration of the relationships among the prediction error matrix E, process variable matrix X, predicted process variable matrix 𝑿̂, and squared prediction error matrix SPE. In Figure 9.7, each row of E contains the row residual or the prediction error for j -th observation (j= 1,2, …J) for all K process variables. Figure 9.8 shows a si milar plot, focusing on the column residual, or the prediction error for each column that represents the k -th process variable (k= 1,2, ...K) in the residual matrix E [13]. Figure 9.8 An illustration of the relationships among the prediction error matri x E, process variable matrix X, predicted process variable matrix 𝑿̂, and column or the prediction error for k -th process variable (column) Each column of E contains the prediction error for one variable. Referring to the discussion of least squares model analysis on pages 165 to 168 of [13], we can find the R2 value fo r the k -th process variable (column) as: 𝑅𝑘2 = 𝑅𝑋,𝑘2=1−𝑉𝑎𝑟 (𝑥𝑘 − 𝒙̂𝒌) 𝑉𝑎𝑟 (𝑥𝑘 ) = 𝑉𝑎𝑟 (𝑒𝑘 ) 𝑉𝑎𝑟 (𝑥𝑘 ) (9.18) The 𝑅𝑘2 value for each process variable will increase with every principal compone nt that is added to the model. The minimum value is 0.0 when there is no principal component and 𝒙̂𝒌=0. The maximum value is 1.0 when we have added the maximum number of principal components with 𝑥𝑘 = 𝒙̂𝒌 and 𝑒𝑘 =0. We can extend th e preceding row residual and column residual concepts to the whole process data matrix X and calculate the R2 value of the entire matrix [13]. This value is the ratio of the variance of X that we can explain with the PCA model over the ratio of variance in itially present in X. R2 = 1−𝑉𝑎𝑟 (𝑿− 𝑿̂) 𝑉𝑎𝑟 (𝑋) =1- 𝑉𝑎𝑟 (𝐸) 𝑉𝑎𝑟 (𝑋) (9.19) By using ML or Python (see Appendix A), or Aspen Technology’s software Aspen ProMV, we can evaluate the R2 value and identify the number of principal components needed to a dequately explain the data variability in X. We have demonstrated this aspect in Appendix A and will illustrate this aspect in our hands -on workshop WS9.1, in which Aspen ProMV shows the R2 value as R2 for different number of principal components. Lastly, page 380 of Dunn [13] explains the concept of determining the number of principal components to use in a model based on cross -validation (CV), originally proposed by Wold [20]. We follow Dunn’s exposition below. The general idea is to divide the process d ata matrix X into 𝐺 groups of rows. These rows should be selected randomly but are often selected in order: row 1 goes in group 1, row 2 goes in group 2, and so on. We can collect the rows belonging to the first group into a new matrix called X(1), and le ave behind all the other rows from all other groups, which we will call group X(−1). So in general, for the 𝑔-th group, we can split matrix X into X(𝑔) and X(−𝑔) . Wold ’s cross -validation procedure asks to build the PCA model on the data in X(−1) using 𝐴 components. Then use data in X(1) as new, testing data. In other words, we preprocess the X(1) rows, calculate their score values, T(1) = X(1)P, calculate their predicted values, X̂︀(1) = T(1)P′ , and their residuals, E(1) = X (1) − X̂︀(1). We repeat th is process, building the model on X(−2) and testing it with X(2), to eventually obtain E(2). After repeating this on 𝐺 groups, we gather up E1, E2, . . . , E 𝐺 and assemble a type of residual matrix, E𝐴,CV, where the 𝐴 represents the number of component s used in each of the 𝐺 PCA models. The CV subscript indicates that this is not the usual error matrix, E. From this, we can calculate a type of 𝑅2 value. We do not call this 𝑅2, but it follows the same definition for an 𝑅2 value. We will call it 𝑄2𝐴 instead, where 𝐴 is the number of components used to fit the 𝐺 models. Q2𝐴 = 1 − Var(E 𝐴, CV)/ Var(X) (9.20) Essentially, Q2 𝐴 is a measure of how well the process variables will be predicted with new data calculated by cross validation. In our han ds-on workshop WS9.1, Aspen ProMV shows the Q2 value as Q2 for different number of principal components with cross validation. 9.1.5 Hotelling’s T2 value from PCA Model In Figure 9.6, we illustrate the score value tk,1 of process variable vector xk on th e first principal component p1 . Let t k,a (k=1,2…K; a=1,2,….A) be the score value of kth process variable xk on the a -th principal component, and s a (a= 1,2,…A) be the variance of the a -th principal component. Then, the Hotelling’s T2 value for the k -th process variable is: T2 = Ʃ (t k,a/sa)2 (9.21) T2 value is a positive, scalar number that summarizes all the score values. It represents the distance from the center of the hyperplane of process variables to the projection of the sample onto the hyp erplane. For samples that are very close to the sample mean gives a T2 value of zero [15]. Figure 9.9 illustrates the concept of Hotelling’s T2 value for an example with two principal components (A=2): T2 = 𝑡12 𝑆12 + 𝑡22 𝑆22 (9.22) Figure 9.9 An illustration of the concept of Hotelling’s T2 value in a two -latent -variable or a two -principal -component space, t 2 versus t 1 In the figure, the equation for T2, Eq. (9.21), is that of an ellipse. T2 expresses how far an observation is from the center of the model in the plane. All points on the ellipse have the same T2 value. We note that