Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Shreeves, Phillip
Master of Science - MSc
University of British Columbia
Raman spectroscopy (RS) is an optical interrogation method used to identify molecules through inelastic light scattering. This process often results in the acquired data having high dimensionality, as each observation has multiple intensities over hundreds of wave numbers. As a result, a common progression in the analysis of the spectra is to first reduce the dimensionality of the data before further examination. This thesis details the use of different dimensionality reduction techniques, including an extension of nonnegative matrix factorization (NMF). NMF is the process of decomposing the data matrix of interest into two lower rank matrices, with the constraint that all matrices must be nonnegative. This algorithm is extended in order to impose constraints on the model by decomposing the one of the specified matrices into two further nonnegative matrices, creating three in total. The application of this technique handles the previously discussed issues while also having added interpretability in comparison to principal component analysis (PCA), which is currently used in RS.
Dimensionality Reduction Techniqueswith Applications in RamanSpectroscopybyPhillip ShreevesB.Sc. Hons., The University of British Columbia, 2018M.Sc., The University of British Columbia, 2020A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe College of Graduate Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)August, 2020c© Phillip Shreeves, 2020The following individuals certify that they have read, and recommend tothe College of Graduate Studies for acceptance, a thesis/dissertation entitled:Dimensionality Reduction Techniqueswith Applications in Raman Spectroscopysubmitted by Phillip Shreeves in partial fulfilment of the requirements ofthe degree of Master of ScienceDr. Jeffrey L. Andrews, I. K. Barber School of Arts & SciencesSupervisorDr. Andrew Jirasek, I. K. Barber School of Arts & SciencesSupervisory Committee MemberDr. W. John Braun, I. K. Barber School of Arts & SciencesSupervisory Committee MemberDr. Sepideh Pakpour, School of EngineeringUniversity ExamineriiAbstractRaman spectroscopy (RS) is an optical interrogation method used toidentify molecules through inelastic light scattering. This process oftenresults in the acquired data having high dimensionality, as each observationhas multiple intensities over hundreds of wave numbers. As a result, a commonprogression in the analysis of the spectra is to first reduce the dimensionalityof the data before further examination. This thesis details the use of differentdimensionality reduction techniques, including an extension of nonnegativematrix factorization (NMF). NMF is the process of decomposing the datamatrix of interest into two lower rank matrices, with the constraint thatall matrices must be nonnegative. This algorithm is extended in order toimpose constraints on the model by decomposing the one of the specifiedmatrices into two further nonnegative matrices, creating three in total. Theapplication of this technique handles the previously discussed issues while alsohaving added interpretability in comparison to principal component analysis(PCA), which is currently used in RS.iiiLay SummaryIn this thesis, a novel algorithm, Group and Basis Restricted NonnegativeMatrix Factorization (GBR-NMF), is developed with motivations from Ramanspectroscopy data sets. Raman spectroscopy is a field in medical physicsthat can identify molecules through an light scattering. These molecules arethen decomposed, using GBR-NMF, into bases that make up said molecules.The GBR-NMF differs from standard NMF due to the fact that suspectedbases and groups can be specified in the model. This model is compared tostandard NMF on simulated data, image data, and the motivating data setsto demonstrate the capability of the GBR-NMF algorithm.ivContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.2 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.3 Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.4 Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . . 3Chapter 2: Background . . . . . . . . . . . . . . . . . . . . . . . 42.1 Principal Component Analysis . . . . . . . . . . . . . . . . . . 42.2 Least Absolute Shrinkage and Selection Operator . . . . . . . 52.3 Logistic Regression . . . . . . . . . . . . . . . . . . . . . . . . 62.3.1 LASSO Logistic Regression . . . . . . . . . . . . . . . 72.4 Nonnegative Matrix Factorization . . . . . . . . . . . . . . . . 72.4.1 Constrained Nonnegative Matrix Factorization . . . . 9Chapter 3: Methodology . . . . . . . . . . . . . . . . . . . . . . 103.1 The Group and Basis Restricted NMF (GBR-NMF) algorithm 103.2 Updating Process . . . . . . . . . . . . . . . . . . . . . . . . . 11vCONTENTS3.3 Proof of Convergence . . . . . . . . . . . . . . . . . . . . . . . 133.4 Pseudocode . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Chapter 4: Applications . . . . . . . . . . . . . . . . . . . . . . . 164.1 Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . 164.2 Facial Expression Data . . . . . . . . . . . . . . . . . . . . . . 184.2.1 Predicting an Unseen Facial Expression . . . . . . . . 204.3 Raman Spectroscopic Cancer Data . . . . . . . . . . . . . . . 214.3.1 Raman Spectroscopic Mouse Data . . . . . . . . . . . 25Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 29Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Appendix A: Derivation of Updating Algorithm . . . . . . . . . . . 35Appendix B: Proof of Algorithm Convergence . . . . . . . . . . . . 38B.1 Proof of Lemma 3.3 . . . . . . . . . . . . . . . . . . . . . . . 38B.2 Proof of Lemma 3.4 . . . . . . . . . . . . . . . . . . . . . . . 39viList of TablesTable 4.1 Summary of the residual sum of squares (RSS) of boththe learned features and scores versus the true featuresand scores for both standard NMF and GBR-NMF. . . 18Table 4.2 Predicting healthy vs. sick using LASSO logistic regres-sion with a 0.5 grade cutoff. . . . . . . . . . . . . . . . 27Table 4.3 Multinomial logistic regression predictions with respectto binning process number 1. . . . . . . . . . . . . . . 28Table 4.4 Multinomial logistic regression predictions with respectto binning process number 2. . . . . . . . . . . . . . . 28viiList of FiguresFigure 2.1 Univariate example of logistic regression with the re-sponse variable being observation class. . . . . . . . . 6Figure 4.1 Row one: Six facial expressions (neutral excluded) av-eraged over the twelve different women. The facialexpressions, in order, are anger, disgust, fear, hap-piness, pain, and surprise. Row two: Six retrievedunconstrained facial features with respect to groupingmatrix corresponding to expressions from the GBR-NMF. Row three: Six retrieved features from the com-pletely unconstrained basic NMF model. . . . . . . . 19Figure 4.2 An original picture of a woman’s face (left) along withreconstructions from the GBR-NMF algorithm (center)and the standard NMF algorithm (right). . . . . . . . 20Figure 4.3 True facial expression (left) in comparison to GBR-NMF prediction of single missing facial expression(center) and GBR-NMF missing facial expression pre-diction using only neutral face (right). . . . . . . . . . 21Figure 4.4 Three averaged Raman spectra displaying three differ-ent cell lines which include lung (red), breast (blue),prostate (green) tumor cells. Grey lines are the ob-served Raman spectra. Unit displayed on the y-axis isarbitrary intensity (AI). . . . . . . . . . . . . . . . . . 22Figure 4.5 Scores on the five (four known/one unknown) chemi-cal bases with respect to the three different cell lines(lung (red), breast (blue), prostate (green)). The Ra-man spectra in the bottom right corner is that of theunknown feature. . . . . . . . . . . . . . . . . . . . . 24Figure 4.6 An experimentally derived Raman spectrum (black)versus a reconstructed spectrum (red) from the analysissummarized graphically in Figure 4.5. . . . . . . . . . 24viiiLIST OF FIGURESFigure 4.7 Example of statistically significant wavenumbers in pre-dicting fibrotic grade. Statistically significant wavenum-bers are shown using both the red and green lines, withred displaying that an increase in intensity results ina decrease in sickness likelihood and green displayingthe opposite. . . . . . . . . . . . . . . . . . . . . . . . 26ixAcknowledgementsMy deepest gratitude goes to Dr. Jeffrey Andrews for taking me on as hisstudent and providing me with guidance over the past 3 years of my academiccareer. I would also like to thank those in the Raman spectroscopy group,led by Dr. Andrew Jirasek, for accepting me in their group and educatingme in their field. I would like to thank Dr. John Braun, Dr. Jirasek, and Dr.Andrews for their contributions on my advisory committee. Finally, I wouldlike to thank both my family and friends for their endless support.xChapter 1IntroductionAs the size of data sets are increasing, the world of data analytics is alsoexpanding. Largely due to the high number of variables, this results in eventhe most commonly known methods to be infeasible at times. This issue hasresulted in an increased focus on a form of data transformation, known asdimensionality reduction, to arise. This is typically done by mapping thedata to a lower dimensional subspace in some way, while also retaining asmuch information as possible from the original data set. There is a largenumber of techniques that are able to perform dimensionality reduction inthis manner, including principal component analysis (PCA, Jolliffe, 2011),vector quantization (VQ, Gray, 1990), factor analysis (Harman, 1960), andnonnegative matrix factorization (NMF, Lee and Seung, 1999). All of thesetechniques follow two key properties of dimensionality reduction specified byWang and Zhang (2012), which includes reducing the dimensionality of thedata and recovering factors that can be interpreted effectively.NMF differs from the other methods specified above primarily due tothe fact that it has the added constraint of nonnegativity throughout itsmatrices. It also benefits as the recovered components are a much more parts-based representation, rather than holistic, which can be beneficial in manyapplications. For example, Lee and Seung (1999) first introduce NMF byapplying it to a data set containing facial images. NMF was able to reduce thedimensionality by identifying components such as the nose, eyes, and mouth.PCA, on the other hand, identified entire facial structures as components,making it difficult to envision the observations as a linear combination of saidparts. The standard NMF algorithm is an unsupervised learning techniqueand is unfortunately sometimes unable to tease out known to exist a priorifactors. We propose a restricted algorithm that allows for one to specifyexpected factors and groupings as added constraints to the model, leading toa more flexible modelling paradigm for the user.The proposed algorithm is particularly useful in the field of Ramanspectroscopy (RS), an optical interrogation method whereby vibrationalmodes of constituent molecules are identified through inelastic light scattering.Raman spectra can provide detailed, multiplexed, information on a range of11.1. Objectivemolecular constituents within a single sample acquisition (Butler et al., 2016;Feng et al., 2017; Pence and Mahadevan-Jansen, 2016). Relative changes inpeak areas and/or positions can be used to identify altered molecular dynamicswithin a system undergoing a given perturbation. Typical dimensionalityreduction practice in the field of RS is to use principal component analysis.However, this is not an ideal model as the constituent underlying spectrashould not be permitted to take on negative values. Furthermore, it isinappropriate to make the assumption that all constituents, which in this caseare represented through the principal components, are uncorrelated. Thisis because the increase of one chemical in a spectrum typically constituteshaving less than another and vice versa. Thus other dimensionality reductiontechniques, such as NMF, are more appropriate for application to thesedata sets. In this thesis, the NMF model is extended to improve overallinterpretation and allow for suspected chemical constituent spectra to bespecified during model fitting. Other models are also applied to said data inorder to further expand knowledge in the field of RS.1.1 ObjectiveRaman spectroscopy techniques typically use principal component analysisin order to reduce the dimensionality of the data. However, these datasets violate some of the assumptions one needs to make. This includes theassumption that components can take on negative values and that componentsare orthogonal of one another. As a result, this is an inappropriate model touse and the data is much better suited for nonnegative matrix factorization.Using NMF adds the constraint that all values must be nonnegative andalso removes the assumption of uncorrelated components. In this thesis, anextension of nonnegative factorization is proposed which allows for previouslyexpected factors as well as known groups to be specified in the model,increasing overall interpretability in comparison to conventional NMF.1.2 Thesis Organization1.2.1 Chapter 2The second chapter reviews modern techniques that will be used andexpanded upon throughout the thesis. Methods such as principal componentanalysis, the least absolute shrinkage and selection operator, logistic regression,and nonnegative matrix factorization are discussed in detail.21.2. Thesis Organization1.2.2 Chapter 3Chapter three expands upon the nonnegative matrix factorization sectiondiscussed in chapter two. It does so by providing a novel extension of thealgorithm through added constraints, allowing for both clustering applicationsas well as the ability to specify suspected factors in the model. These imposedconstraints, along with the algorithm’s updating process and optimizationcriterion, are outlined in this section.1.2.3 Chapter 4This chapter provides proof of concept of the algorithm provided in chapterthree. It shows applications of the algorithm on multiple data sets and alsodisplays other dimensionality reduction techniques that can be applied inRaman spectroscopy examples.1.2.4 Chapter 5The final chapter of the thesis concludes and summarizes the work thatwas completed throughout. The impact that this thesis provides and futurework to be performed are also discussed.3Chapter 2BackgroundHerein, several common and modern techniques for dimensionality reduc-tion and modelling are reviewed.2.1 Principal Component AnalysisOne of the most commonly used dimensionality reduction techniques,Principal Component Analysis (PCA) was first developed by Pearson (1901)and is used in many different situations and some primary goals of the modelinclude variable selection, classification, and dimensionality reduction (Woldet al., 1987). PCA uses linear combinations of the original variables in thedata set to generate new variables. The new variables created are uncorrelatedand are made to retain as much variance from the original data as possible,while also reducing dimensionality (Jolliffe and Cadima, 2016). This bycalculating the first principal component (PC) by achieving the maximumpossible variance, then calculating the second PC in the same fashion whileensuring it is uncorrelated with the first (Jolliffe, 1990). The process ofcalculating components continues until the number of PCs matches the totalnumber of variables in the original data set. In particular, the weights ofthe linear combination for each of the new variables (PCs) are found via theeigenvectors of the covariance/correlation matrix of the original variables.From there, the number of PCs to be kept can be determined in regards tohow much variance one desires to keep. The eigenvalues of said eigenvectorsare then proportional to the percentage of variance from the original dataassociated with that PC. Each observation is then mapped to the PC spacevia weights on each original variable. Notably both the eigenvectors andweights can take on negative values. While removing correlated features isnot desirable in certain applications, it often is a benefit of PCA. It is alsoadvantageous to be able to choose how much information one wants to retainwith respect to choosing a number of principal components. The caveatregarding choosing said number is that information is guaranteed to be lostunless all components are retained. Along with loss of information, PCA alsobecomes difficult to interpret compared to the original variables in complex42.2. Least Absolute Shrinkage and Selection Operatorapplications. For more comprehensive information on Principal ComponentAnalysis, see Jolliffe (1986), Abdi and Williams (2010), and Jolliffe (2011).2.2 Least Absolute Shrinkage and SelectionOperatorThe “Least Absolute Shrinkage and Selection Operator”(Tibshirani, 1996),more commonly referred to as LASSO, is a well known technique used forvariable selection and shrinkage of coefficients in generalized linear models.First developed with respect to the commonly known ordinary least-squareslinear regression model, this method has expanded to other familiar modelswhich includes the Cox regression (Tibshirani, 1997) and the logistic regression(Shevade and Keerthi, 2003). In an application of the ordinary least squaresregression, coefficients are estimated by minimizing the following function(α,β) = argminα,βn∑i=1(yi − α−p∑j=1(βjxij))2. (2.1)However, commonly known issues with said algorithm include the largeamounts of variability between estimated variables and true values and lackof interpretability with respect to high dimensionality models (Tibshirani,1996). These issues can be resolved through the addition of a penalty term,which is added to the objective function as shown below(α,β) = argminα,βn∑i=1(yi − α−p∑j=1(βjxij))2 + λ || β ||1 (2.2)where α and β represent the intercept and model coefficients respectively.The additional term is effectively equivalent to minimizing Equation 2.1using a constraint (∑ | βj |< s) on the size of the coefficients in the model(Tibshirani, 2011). This, in turn, allows for the adjustment of said λ termin order to control the total amount of shrinkage of the overall model. Forexample, as λ approaches infinity, it is notable that the beta terms would needto be shrunken to zero in order to optimize the function. Contrarily, if λ takeson a very small value, then the coefficients are not substantially restricted inmagnitude which may result in large values. Additional information regardingLASSO can be found in sources by Kwon et al. (2013) and Hastie et al. (2015).52.3. Logistic RegressionFigure 2.1: Univariate example of logistic regression with the response variablebeing observation class.2.3 Logistic RegressionLogistic regression is a classification technique used to predict an observa-tion’s class with the variables at hand and determine the relationship betweenthem. The model in a two group scenario is characterized using the followingequationP (Y = 1 | xi) = 11 + e−(α+x′iβ)(2.3)where P (Y = 1 | xi) represents the probability of the observation belongingto the class of interest. This causes said model to be frequently used in thefield of epidemiology due to P (Y = 1 | xi) taking a probabilistic value. Inunivariate cases, this model is well known to have an S-shaped probabilitycurve as shown in Figure 2.1 which allows for a distinct cut-off point betweenclasses. This technique can further be expanded to a multi-class scenario(Böhning, 1992) that takes similar form to its binomial counterpart. Theprobability of observation i belonging to group g ∈ {1, ..., G} can be calculatedin the subsequent mannerP (Y = g | xi) = eα+x′iβ(g)1 + eα+x′iβ(g)for g = 1, ..., G− 1 (2.4)P (Y = G | xi) = 11 + eα+x′iβ(G)(2.5)62.4. Nonnegative Matrix Factorizationwhere β(g)j represents the jth coefficient in the gth class model and the Gthgroup represents the reference variable. Further information on this topic canbe found in sources by Kleinbaum et al. (2002), Hilbe (2009), and Hosmer Jret al. (2013).2.3.1 LASSO Logistic RegressionUnlike the conventional LASSO regression, LASSO logistic regressionimposes a constraint on the negative log-likelihood of the model. Thislikelihood is given by Huang et al. (2008) as`n(β) = −n∑i=1(Yi log pi(x′iβ) + (1− Yi) log[1− pi(x′iβ)]) (2.6)where pi(t) = et/(1 + et) and is constrained such thatβˆ(λ) = argmin(n−1`n(β) + λ || β ||1) (2.7)by Fonti and Belitser (2017). This constraint is added in a similar fashionwhen LASSO is applied to the multinomial logistic regression (Zhai, 2018).2.4 Nonnegative Matrix FactorizationNonnegative matrix factorization was originally developed by Lee andSeung (1999) as a method of finding the major underlying factors of a modeland seeing how each individual observation is composed of the factors. Thefactors are found by mapping the observations to a latent factor space thataccounts for the vast majority of variance in the model. Each observation isthen written as a linear combination of the discovered factors. These factorsand coefficients of the linear combinations are found by decomposing thenonnegative data of interest X ∈ Rn×p≥0 (where n is the number of observationsand p is the number of variables) into two lower rank nonnegative matricesW ∈ Rn×q≥0 and H ∈ Rq×p≥0 such thatX ≈WH (2.8)where with q underlying factors, we can describe the span of the columns ofH as the q-dimensional subspace and W as the scores for each observation onthat subspace. Thus, each observation xi can be viewed as an additive linearcombination of the features contained in the rows of H that are weighted by72.4. Nonnegative Matrix Factorizationthe scores contained in the columns of W. This equation is an approximationof the assumed true decomposition of XX =WH+ ε (2.9)with ε ∈ Rn×p represents the residual or noise matrix of the model. This noiseε = X−WH, is a key component in the model as its assumed probabilitydistribution assists in determining what objective function the data needsto be optimized with respect to (Wang and Zhang, 2012). For example,assuming a Gaussian distributed residue matrix, the Frobenius normDF (X |WH) = 12|| X−WH ||2F=12∑ij(xij − [WH]ij)2 (2.10)is the function that is to be optimized (Lee and Seung, 2001; Wang and Zhang,2012). Another function commonly used is the generalized Kullback-Leiblerdivergence (GKLD) which is used when assuming Poisson distributed noise.For this thesis, we will assume Gaussian noise as the simulated data sets werecreated with said noise and preliminary diagnostics regarding real data setsfor complex Raman spectra such as those in Sections 4.3 and 4.3.1 did notreject the assumption of Gaussian noise.Lee and Seung (2001) describe that NMF can be formulated as twooptimization problems, which are to minimize both || X −WH || andD(X | WH) simultaneously and while these two functions are convex inW or H individually, they are unfortunately not convex together. Becauseof this, finding global minima is not guaranteed and algorithms which findlocal minima are predominantly considered at the moment. The prototypicalalgorithm for finding these matrices is an alternating multiplicative updatingprocess that uses the following updates:hij ← hij (WTX)ij(WTWH)ij, wij ← wij (XHT )ij(WHHT )ij. (2.11)This updating process is considered slow by most accounts and scales dra-matically in elapsed time when requiring both the rank and dimension of Xto increase. Suggestions to increase the performance of the algorithm includeapplying gradient descent algorithms (Lin, 2007) or use of the conjugategradient (Zdunek and Cichocki, 2007). While NMF is advantageous dueto it’s added interpretability in comparison to PCA, it does fall victim toslower computation speed (Devarajan, 2008). Aside from the assumption ofnonnegativity, Lee and Seung (1999) also note that no further assumptionsare made about the statistical independencies of the underlying factors.82.4. Nonnegative Matrix Factorization2.4.1 Constrained Nonnegative Matrix FactorizationVarious authors have considered constraints to address questions concern-ing clustering and classification with the NMF algorithm. Ding et al. (2005)suggest that when the Frobenius norm is used as the objective function, thestandard NMF algorithm bears resemblance to a relaxed form of K-meansclustering (Hartigan and Wong, 1979) where one factor matrix contains thecentroids of the clusters and the other contains indicators of cluster member-ships. This is then expanded upon by Li et al. (2007) in a semi-supervisedclustering manner, where limited knowledge of cluster membership is present.Wagstaff et al. (2001) specify how semi-supervised clustering is an algorithmperformed with respect to two separate constraints. These two constraintsare that there is the existence of must-link observations such that said ob-servations must be grouped into the same cluster and that theres also existscannot-link observations that are required to be in separate clusters.When this knowledge is limited to a number of observations l < n, asemi-supervised clustering algorithm is needed. Liu et al. (2011) addressthese constraints by introducing the decomposition of the feature scoresmatrix W into two lower rank matrices, consisting of a constrained groupingmatrix as well as an ‘auxiliary matrix’. Here, the grouping matrix uses thea priori label information and stores it in the form of an indicator matrixC = [ c1 c2 ... cg ] where g is the total number of classes. This matrix is thenexpanded upon with each unsupervised observation being assigned its ownclass such that W =[C 00 I]. Both W and the auxiliary matrix are thenupdated iteratively with that information taken into account. While thisalgorithm serves a different purpose than the one that will be described inChapter 3, it still is able to provide a contribution to the updating processdescribed in Section 3.2. Yang et al. (2018) further expand upon the workof Liu et al. (2011) by adding in a sparsity constraint, making it a dualconstrained algorithm.9Chapter 3Methodology3.1 The Group and Basis Restricted NMF(GBR-NMF) algorithmThe conventional NMF algorithm proposed by Lee and Seung (1999) isa completely unsupervised algorithm which is cause for concern in certainapplications sometimes. This is because one cannot specify any suspectedfactors in the model or specify any clustering constraints if there are suspectedgroups within the data. Here we propose a variation of the work by Liu et al.(2011) by similarly utilizing a further decomposition of the data of interestinto matrices W, A, and S such thatX =WAS+ , (3.1)where it is apparent that H from Equation 2.8 is further decomposed into aq × q auxiliary matrix A and the q × p matrix S containing the underlyingfactors. However, this algorithm differs in the way that restrictions areimposed on the matrices. Liu et al. (2011) restrict their W matrix as statedin Section 2.4.1, whereas the new algorithm restricts a number of columns inW through the use of a complete indicator matrix and also restricts somerows in S. The added restrictions allow for one to specify known groups in Wand suspected factors in S, further improving the interpretability of the model.It is important to note that while there is further decomposition, the matricesW and S can still be interpreted the same way as W and H respectively.However, due to the fact that we will now be imposing a constraint on S, A isneeded to perform adjustments on the constrained factors. We propose thatA be initialized as an identity matrix which, through the iterative updatingprocess, will act as a scaling matrix to the rows of S — a detail which willbe shown to be quite significant in Section 4.3.As previously specified, constraints will now be added to S, herein creatinga restricted model. Assuming that there are q underlying factors in the model,suppose there exist k known and q − k unknown factors in the model. Thenthese k known factors can be specified as factors in S that are not updated103.2. Updating Processthrough the iterative process, while the q − k factors are still included in theupdating process. This allows for the algorithm to take into considerationthe factors that a user may know to be present and use them to estimateother factors that may not be known to the users.Branching from the work of Ding et al. (2005), an optional constraintregarding the W matrix has also been added for applications regardingclustering. As previously described in Section 2.4.1, basic NMF resemblesa soft K-means clustering algorithm under use of the Frobenius norm as anobjective function. Here, the features in S represent the centroids in thealgorithm and the scores in W are how each observation scores with respectto said centroids. With this in mind, we have allowed the constraint of thefirst g columns of W to be a known grouping matrix. These columns arenever updated in the model and the factors regarding these columns are keptunconstrained.Thus, in the case where both constraints are employed we have g + k ≤ qwith the first g columns of W and rows of S responsible for the clusteringresults with the next k accounting for the constrained factors. This allows forthe major differences in groups to be found while the minor dissimilaritiesare found in the other q − g factors.3.2 Updating ProcessAssuming we are updating with respect to Gaussian distributed error, weare required to update with respect to the Frobenius normDF (X |WAS) = 12|| X−WAS ||2 (3.2)∝ || X−WAS ||2 (3.3)as the objective function of interest with the constraint that wij ≥ 0, aij ≥ 0,and sij ≥ 0 for all ij ∈ R. As previously stated in Section 2.4, it is notfeasible to find global minima as the GBR-NMF is not convex in W,A, andS, which requires an updating process to instead find local minima. To ensurea suitable local minimum has been achieved, the GBR-NMF algorithm wasapplied to different data sets multiple times in order to ensure consistency inthe accuracy of the algorithm. These different test cases included differentinitializations of the W, A, and S matrices in an attempt to find differentlocal minima. If one finds themselves in the situation where they believe thelocal minimum is unsatisfactory, they too could use the algorithm in multiplecases and compare Frobenius norms to find the most suitable answer.113.2. Updating ProcessApplying the trace function, as well as both the cyclic and transposeproperties associated with it, we haveD(X |WAS) ∝ Tr((X−WAS)(X−WAS)T )= Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT ).(3.4)However, this does not account for the nonnegativity constraints previouslyspecified. In order to do so, a Lagrange function (L) with multipliers αij , βij ,and γij are required for constraints wij ≥ 0, aij ≥ 0 and sij ≥ 0 respectively(note that α = [αij ], β = [βij ], and γ = [γij ]). This gives the followingLagrange functionL =Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )+ Tr(αWT ) + Tr(βAT ) + Tr(γST ).(3.5)Now needing the derivatives of L to be zero with respect to W, A, and S wehave∂L∂W= −2XSTAT + 2WASSTAT + α = 0,∂L∂A= −2WTXST + 2WTWASST + β = 0, (3.6)∂L∂S= −2ATWTX+ 2ATWTWAS+ γ = 0.Then applying the Karush-Kuhn-Tucker conditions αijwij = 0, βijaij = 0,and γijsij = 0, we obtain(XSTAT )ijwij − (WASSTAT )ijwij = 0,(WTXST )ijaij − (WTWASST )ijaij = 0, (3.7)(ATWTX)ijsij − (ATWTWAS)ijsij = 0.Which leads to the following updating ruleswij ← wij (XSTAT )ij(WASSTAT )ij, (3.8)aij ← aij (WTXST )ij(WTWASST )ij, (3.9)sij ← sij (ATWTX)ij(ATWTWAS)ij. (3.10)Further detail of the derivations in this section can be found in AppendixA.123.3. Proof of Convergence3.3 Proof of ConvergenceIn order to prove that the above updates are correct, the objective functionmust be non-increasing according to the updates. The objective function isthen invariant under the updates if and only if W, A, and S are at stationarypoints (Liu et al., 2011). In order to prove this, we must borrow a lemmastemming from updates to the EM algorithm (Dempster et al., 1977; Wu,1983) provided below.Lemma 3.1. If an auxiliary function, G, exists for F (x) and satisfies theconditions1. G(x, x′) ≥ F (x)2. G(x, x) = F (x)then F is non-increasing using the updatext+1 = argminxG(x, x′)Given that Fwij is defined as the parts of equation 3.2 relevant to wijand F ′wij is the derivative of such, we then need to prove the following threelemmas.Lemma 3.2. The functionG(w,wij) = Fwij (wtij) + F′wij (wij)(w − wtij) +(WASSTAT )ijwtij(w − wtij)2is an auxiliary function to Fwij .Lemma 3.3. The functionG(a, aij) = Faij (atij) + F′aij (aij)(a− atij) +(WTWASST )ijatij(a− atij)2is an auxiliary function to Faij .Lemma 3.4. The functionG(s, sij) = Fsij (stij) + F′sij (sij)(s− stij) +(ATWTWAS)ijstij(s− stij)2is an auxiliary function to Fsij .133.4. PseudocodeThese lemmas are proved in a very similar fashion. As a result, lemma 3.2will be proved below, while Lemmas 3.3 and 3.4 can be found in appendix B.With minor modifications of the proof by Liu et al. (2011), the following canbe proved.Proof. It is clear that G(w,w) = Fwij (w). According to the definition of anauxiliary function, it only needs to be shown that G(w,wtij) ≥ Fwij (w). Thiscan be done using the Taylor series expansion of Fwij :Fwij (w) = Fwij (wtij) + F′wij (w − wtij) +12F ′′wij (w − wtij)2with F ′′wij being the second order derivative of Fwij . We can show thatFwij = (∂D∂W)ij = (−2XSTAT + 2WASSTAT )ij ,F ′′wij = (∂2D∂W2)ij = 2(ASSTAT )jj .As F ′′wij has now been solved for, it is equivalent to prove(WASSTAT )ijwtij≥ 12F ′′ij = (ASSTAT )jj (3.11)in order to reach Lemma 3.2. Now we have(WASSTAT )ij =q∑`=1(W)i`(ASSTAT )`j (3.12)≥ wtij(ASSTA)jj (3.13)which when rearranged proves Equation 3.11.3.4 PseudocodeFor the observed n× p data matrix, X, with k suspected factors in S andq total factors, the algorithm can be summarized using the following steps:1. If an n×c cluster membership matrix is specified, set the first c columnsof the matrix, W, to be said matrix. Initialize the remaining n× (q− c)elements sampling from the uniform distribution U(min(X),max(X)).Otherwise, initialize the entire matrix by sampling from the uniformdistribution specified above.143.4. Pseudocode2. Initialize A such that it takes the form of a q × q identity matrix.3. Initialize the remaining (q−k) columns of S sampling from the uniformdistribution on the bounds of the known factor values.4. Calculate D(0)F using equation 3.2.5. Set z = 1 and enter the updating process:(a) Update the unspecified elements of W using Equation 3.8.(b) Update all elements of A using Equation 3.9.(c) Update the unspecified elements of S using Equation 3.10.(d) Calculate D(i)F using equation 3.2.(e) If D(z−1)F −D(z)F is less than a specified δ value or the value of i isequal to the maximum iteration count, move to step 6. Otherwiseset z = z + 1 and return to step 5a.6. Return the updated values of W, A, and S.It is important to note that stopping criterion of step 5e can be determinedin a variety of ways. GBR-NMF provides two possible methods which includelack of progress and a maximum iteration count. Unless otherwise specified,the maximum iteration count is used as lack of progress depends greatly onthe value of δ provided. Data sets with large variability may not be able toattain a small epsilon value. Contrarily, if the specified δ value is too largethen the stopping criterion will be achieved too quickly, leaving room forimprovement. To battle this, a high maximum iteration count of 50,000 isset to ensure a local maximum is reached.15Chapter 4Applications4.1 Simulated DataA simulation was generated in order to test the accuracy of the algorithmusing n = 400 observations, p = 2000 variables, g = 4 groups, and q = 7true underlying factors. This simulation was created to test the accuracywith respect to the clustering aspect of the algorithm as well as how it doespredicting scores on known factors. The simulation stipulates that each groupdifferentiates itself from the others by having one unique factor that theothers did not. Each one of the four groups had a unique factor, along withthree other factors that were common with the other groups. Thus, whenrunning the algorithm with respect to this data, the first four columns ofthe W matrix are set to be the known classification matrix of the model. Toadd to the simulations, one of the three common features in S were held as aconstrained factor to determine how well it’s true scores were being found.The matrices below further illustrate which columns and rows in the modelare being held constant (in bold) throughout the updating process.W =1 0 0 0 w1,5 w1,6 w1,70 1 0 0 w2,5 w2,6 w2,70 0 1 0 w3,5 w3,6 w3,70 0 0 1 w4,5 w4,6 w4,7.1 0 0 0 w5,5 w5,6 w5,7.....................0 0 0 1 w400,5 w400,6 w400,7(4.1)A =a1,1 0 0 0 0 0 00 a2,2 0 0 0 0 00 0 a3,3 0 0 0 00 0 0 a4,4 0 0 00 0 0 0 a5,5 0 00 0 0 0 0 a6,6 00 0 0 0 0 0 a7,7(4.2)164.1. Simulated DataS =s1,1 s1,2 . . . s1,2000s2,1 s2,2 . . . s2,2000s3,1 s3,2 . . . s3,2000s4,1 s4,2 . . . s4,2000s5,1 s5,2 . . . s5,2000s6,1 s6,2 . . . s6,2000s7,1 s7,2 . . . s7,2000(4.3)Entries for W, A, and S were all created using random uniform distributionsand then scaled accordingly such that the column sums of W are equalto n and the area under the curve of the columns of S were equal to 1.Thus the diagonal elements of A were scaled inversely such that the overallreconstructions would remain the same as those prior to scaling.These data sets were created and decomposed via the GBR-NMF 100different times in order to ensure accuracy of the results. The residual sumof squares (RSS) comparing the true S with the estimated S were thencalculated and a mean and standard error of these numbers were takenacross all features (found in Table 4.1). This was similarly done with theresulting scores matrix W and auxiliary matrix A to quantify the algorithm’soverall estimation performance. More specifically, to determine how close thevalues of the calculated matrices are to the true values. Table 4.1 displaysthat the algorithm does exceptionally well in predicting the features in bothsimulations. Meanwhile, the accuracy of the scores matrix was also quiteaccurate in comparison to that of the RSS regarding the A matrix. This isprimarily due to the scalings on the two outside matrices which causes thevalues to become quite close. Therefore, a more fitting comparison would beto the unconstrained model, which is also displayed in Table 4.1. This showsa drastic increase in both the accuracy of the scores and factors in comparisonto the original model. It is important to note that this is a direct comparisonbetween the standard NMF model and GBR-NMF. Comparisons to PCA orthe semi-supervised NMF model of Liu et al. (2011) would be inaccurate asa direct comparison is not available. Principal components acquired fromPCA are uncorrelated with respect to one another, making it difficult tomake a direct comparison to the true factors created. Similarly, the differentconstraints in the work of Liu et al. (2011) create different interpretations ofthe factor matrices.174.2. Facial Expression DataW constraint held Standard NMFAverage Standard ErrorW RSS 4,786,092,873 19,172,015,537H RSS 0.000040 0.000010W constraint held GBR-NMFAverage Standard ErrorWA RSS 2,613,972,969 731,945,188S RSS 0.000008 0.000004Table 4.1: Summary of the residual sum of squares (RSS) of both the learnedfeatures and scores versus the true features and scores for both standardNMF and GBR-NMF.4.2 Facial Expression DataA common research application of nonnegative matrix factorization isdata regarding facial expressions (Buciu and Pitas, 2004; Guillamet andVitria, 2002; Zhi et al., 2010) typically for the identification of certain facialexpressions or people. The pain expression data set (Psychological ImageCollection at Sterling, 2019) consists of 84 observations of 241× 181 pixelfacial images. These observations can be further subsetted into the faces of12 different women; each of which displaying 7 different facial expressions.An average face of each facial expression (aside from neutral) is displayed inrow one of Figure 4.1. These six facial expressions are the grouping variableof interest in this data set and are placed as a grouping constraint in theW matrix in the GBR-NMF. Along with said grouping constraints, twelvefeatures were also constrained to be the neutral facial expressions of thedifferent women. These constraints are determined to test what bases areacquired from the grouping variables as well as how close one can get to thetrue facial expression using a person’s neutral face and grouping information.From these neutral faces, the other bases are recovered and built on top of thisneutral face in order to recover the original images. The six bases recoveredusing GBR-NMF resemble key facial features explanatory of each emotion.Row two of Figure 4.1 identifies the unconstrained features retrieved fromthe model. The fourth image of the second row displays a smile expression,which is portrayed well by highlighting the teeth and perked cheeks. Theseresults are more interpretable and separated than those that are recoveredfrom the standard NMF algorithm as this algorithm tends to recover featuresthat appear to more closely resemble particular women’s faces rather thanemotion-specific representations. The standard NMF results are shown in184.2. Facial Expression DataFigure 4.1: Row one: Six facial expressions (neutral excluded) averaged overthe twelve different women. The facial expressions, in order, are anger, disgust,fear, happiness, pain, and surprise. Row two: Six retrieved unconstrainedfacial features with respect to grouping matrix corresponding to expressionsfrom the GBR-NMF. Row three: Six retrieved features from the completelyunconstrained basic NMF model.194.2. Facial Expression DataFigure 4.2: An original picture of a woman’s face (left) along with recon-structions from the GBR-NMF algorithm (center) and the standard NMFalgorithm (right).row three of Figure 4.1 which displays the six features of an unconstrainedNMF sorted to best align with the original expressions.This application displays that we are acquiring better separation in theGBR-NMF in comparison to the basic NMF with respect to the recoveredfacial expressions. These interpretability benefits could be assumed to resultin a decrease in overall model performance. However, as demonstrated inFigure 4.2, the reconstructions still do compare quite well to those of theunconstrained NMF model, making this a worthy trade-off.4.2.1 Predicting an Unseen Facial ExpressionIn this section, we consider the possibility of generating previously unseenimages using the GBR-NMF algorithm. Suppose one wanted to predict aperson’s facial expression that is not included in the original model fit. Thenthe GBR-NMF could be used in order to predict what that expression mightlook like. The only information that would be needed in order to do so arethe acquired bases from the model, as well as one’s neutral facial expression(ie. the S matrix). This problem was tackled in two different ways. Thefirst scenario involves simply removing one observation in the data set andattempting to predict the face that was removed. The second requires allof the data regarding a person’s facial expressions to be removed from themodel and then predicting what their expression would look like; again, usingonly the neutral facial expression and the bases acquired. The predictionsof these faces were made by scoring 1 on the desired facial expression and1 on the neutral face, and using this linear combination as the prediction.Acquisitions of the predicted faces, along with the true facial expression are204.3. Raman Spectroscopic Cancer DataFigure 4.3: True facial expression (left) in comparison to GBR-NMF predic-tion of single missing facial expression (center) and GBR-NMF missing facialexpression prediction using only neutral face (right).shown in Figure 4.3. While it is clear that these reconstructions are not asclose to the original image as the results in Figure 4.2, the added benefit ofthe model is interpretability. In comparison to the conventional NMF, a newfacial expression not seen by the algorithm can be predicted due to how easythe factors are to understand. The conventional NMF algorithm would notbe able to to this in many cases due because it’s factors cannot be interpretedas easily.4.3 Raman Spectroscopic Cancer DataIn this example, Raman spectroscopy can be used to provide detailedinformation on the changes in relative concentrations of classes of molecularcompounds (proteins, lipids, DNA, metabolites, etc.) within cellular andtissue environments exposed to ionizing radiation used in cancer therapy(Matthews et al., 2015; Harder et al., 2015; Matthews et al., 2010; Paidiet al., 2019). Matthews et al. (2015) introduce data containing 3240 spectracontaining 582 Raman intensities (arbitrary units) at different wavenumbers(cm−1). These spectra are separated into three different groups based oncell type, specifically lung (H460), breast (MCF-7), and prostate (LNCaP)tumours. The spectra, overlaid with an average of each cell type, can beseen in the Figure 4.4. In Raman spectroscopy, chemical constituents areidentified through comparison with known chemical constituent spectra andliterature line lists exist for the most common cellular components. Matthewset al. (2015) performed a principal component analysis on the spectra andfound that the first component closely resembled the metabolite glycogen, a214.3. Raman Spectroscopic Cancer DataFigure 4.4: Three averaged Raman spectra displaying three different celllines which include lung (red), breast (blue), prostate (green) tumor cells.Grey lines are the observed Raman spectra. Unit displayed on the y-axis isarbitrary intensity (AI).key component that is produced in radioresistant cells’ response to radiationdamage. The aim of this study is to further expand on the work of Matthewset al. (2015) by recovering more biological basis using GBR-NMF. Using theGBR-NMF allows for the constraint of the factor matrix, allowing for thepreviously recovered glycogen to be specified in the model. Further recoveredbases are then analyzed to discover new biological information with respectto the chemical composition of the cells in the data set.There are two main, and related, disadvantages to using PCA in thiscontext. Firstly, permitting negative scores and component coefficientsviolates the physical understanding of the chemical composition of a cell —reconstructions of the spectra could be built in a subtractive manner wherean overabundance of one chemical is estimated and then subtracted viaa negative score on another chemical. Secondly, and perhaps related tothe problems with previous assumptions, Matthews et al. (2015) found alack interpretability in additional principal components of the model. Thestandard NMF algorithm gains a key advantage over PCA due to the factthat it is forced to be an additive model with respect to chemical constituents,which matches the physical construction of Raman spectroscopic data. Denget al. (2020) further illustrates this point by performing NMF on the same224.3. Raman Spectroscopic Cancer Datadata set and not only discovering glycogen, but also a spectrum very similarto that of a lipid. While using standard NMF is beneficial over PCA, useof the GBR-NMF allows making use of previous results from Matthewset al. (2015) wherein an important chemical, glycogen, was discovered andbiologically confirmed. As such, we utilize GBR-NMF by constraining anacquired spectra of pure glycogen within S, and further constrain W withthe known cell lines in the data. We permit one further basis row and groupcolumn to be fully unconstrained.Figure 4.5 provides a scatterplot matrix of the scores from W on the offdiagonal and the recovered spectral bases along the diagonal. The first threebases along the diagonal were recovered under group constraints on the threetumour types. We note that all three bases in these cases closely approximatea standard cell spectra, with slight variations among some wavelengths. Theglycogen scores for the green group (prostate, LNCaP) fall below a score of 1.5;which matches the result put forward by Matthews et al. (2015) that LNCaPprostate tumours exhibit lower glycogen production in response to radiationthan the other two cell lines (H460, MCF7). Furthermore, LNCaP prostatetumours are also known to be more sensitive to radiotherapy than thoseof H460 and MCF-7. Interestingly, the estimated results from GBR-NMFshow that the prostate cells are not fully differentiable from the other celllines based on glycogen alone, as there is substantial overlap in the glycogenscores, even if the LNCaP cells have lower, on average, scores. However, if weconsider the fifth, fully unconstrained basis that was recovered (fifth diagonalalong Figure 4.5), we also see lower scores on average from the LNCaPcells. In fact, if the scores are viewed in tandem between glycogen and thiscurrently unknown spectral component, we see complete separation betweenLNCaP cells (green) and the H460 and MCF-7 cell lines (blue and red). Thus,applying the GBR-NMF method has provided a new spectral componentto investigate with regards to potential chemicals that could be indicativeof cell radiosensitivity. The current result presented here allows for a moredirect focus on individual biochemical cellular response to radiotherapy ascompared with the original basic PCA analysis. Further work is required toelucidate the biochemistry of radiation response and is beyond the scope ofthis thesis.A natural question that arises is how well the model estimated by GBR-NMF approximates the original spectra. Consider the plot in Figure 4.6showing two spectra, one being a randomly selected spectrum from theoriginal data and the other being a reconstruction of that spectrum fromGBR-NMF. From a spectroscopic viewpoint, this would be considered anexcellent reconstruction.234.3. Raman Spectroscopic Cancer DataFigure 4.5: Scores on the five (four known/one unknown) chemical bases withrespect to the three different cell lines (lung (red), breast (blue), prostate(green)). The Raman spectra in the bottom right corner is that of theunknown feature.Figure 4.6: An experimentally derived Raman spectrum (black) versus areconstructed spectrum (red) from the analysis summarized graphically inFigure 4.5.244.3. Raman Spectroscopic Cancer Data4.3.1 Raman Spectroscopic Mouse DataThe data set in this example contains 4074 spectra containing 672 differentRaman intensities. These spectra were acquired from slices of lungs ofmice that had varying degrees of fibrosis and pneumonitis. The degree ofsickness in these observations were assessed in accordance with previouslyestablished methodologies and were recorded. It is important to note thateach observation does not hold its own unique sickness gradings, but insteadgradings of the region that the spectrum lies in. These gradings are theproportion of the region that are considered sick and each observation takesthese values due to the fact that it is not possible to locate exactly where thespectrum was acquired in the region of interest.Of the 4074 observations obtained, not all are of the lung tissue describedabove. This is because some of these spectra taken hit air pockets in thelungs, instead returning the media in which the tissue sample is encased.In order to filter the data set such that only desirable spectra are retained,the GBR-NMF was applied to the data using three different bases - ‘goodspectra’, ‘bad spectra’, and media. Spectra considered ‘good’ are thosewhich are representative of Raman spectra if tissue, whereas the media isthe solution which the tissues were immersed in during the acquiring of thespectra. The ’bad spectra’ are those which resembled the cells to a certainextent, but still displayed peaks which were present in the media. The baseswere created by averaging 10 observations considered to take that category;for example, to create the ‘good’ basis, 10 desirable observations were foundin the data and averaged to create said basis. This was also done for boththe ‘bad’ basis, while the media was a pure spectrum of the chemical. Uponinitial analysis, a distinct cutoff was determined where any observations witha media score greater than 1 or a ‘good’ spectrum score less than 1 is removedfrom the data set. As a result, 2415 of the original 4074 observations werekept. It is important to note that this is not considered outlier removal as notonly observations which differ significantly from the rest are being removed.Instead, one should consider it to be group removal, as all spectra which donot fit the desired chemical composition are being removed.The primary goal with this data was to be able to predict the degreeof sickness in a lung tissue sample. More specifically, if one were given aRaman spectrum, could the cell corresponding to that observed spectrumbe predicted as sick or not. This can be done through the use of a LASSOlogistic regression model, which uses the observation intensities to predictsickness. Based on this model, certain wave numbers are chosen as significant,while others’ coefficients are shrunken to zero in order to simplify the model254.3. Raman Spectroscopic Cancer DataFigure 4.7: Example of statistically significant wavenumbers in predictingfibrotic grade. Statistically significant wavenumbers are shown using both thered and green lines, with red displaying that an increase in intensity resultsin a decrease in sickness likelihood and green displaying the opposite.as a whole (see figure 4.7). However, “sickness” is not an actual variable inthe data at hand. Each region of observations taken instead has a specificgrade between zero and one given to it, which represents the proportion ofthe observation considered infected. As a result, a threshold needs to bechosen in order to create the response variable of the model. The thresholdchosen was 0.5, meaning all observations in a region with a grade below 0.5is considered healthy and those in regions with grades above are consideredsick. This was the initial cutoff point chosen as it was considered the mostrational value when choosing a threshold. If more than 50% of the region isconsidered infected, then all cells are considered disease, causing separationof regions according to the majority of the regions’ cell health. While thismethod is shown to be effective, it is not quite the correct model to be using.Suppose a region has a grading of 0.6. This would then mean that roughly40% of the observations sampled from this region would not have the disease.Therefore, to classify all observations in this region as infected would befalse. Instead, a classification process that has classes other than zero andone should be considered. Herein, multinomial logistic regression is used toclassify observations in a multi-class scenario very similarly to the process264.3. Raman Spectroscopic Cancer DataTable 4.2: Predicting healthy vs. sick using LASSO logistic regression with a0.5 grade cutoff.PredictedTrue Healthy SickHealthy: 1654 7Sick: 4 750discussed before. This case, however, considers classifying according to theactual grades provided as opposed to sick/healthy. In turn, this then providesa more accurate diagnostic of the observations’ level of sickness. It should benoted that this method does not provide an absolute prediction as shown withthe previous logistic regression model. For example, suppose an observationis predicted to have a grade of 0.3. Then the accurate interpretation wouldbe that it is predicted that 30% of the cell’s region is predicted to be sick.This then could be further interpreted as the cell having 30% probability ofbeing infected. As the grading process can be somewhat subjective, it shouldalso be noted that these grades could be binned into different risk levels. Thegrades at hand (0,0.05,0.1,0.3,0.5,0.6,0.9,1) can be separated using the twodifferent binning methods outlined below1. Low (0-0.1), Medium (0.3-0.6), High (0.9-1)2. Zero (0), Low (0.1-0.3), Medium (0.5-0.6), High (0.9-1)These methods allow for a more regulated approach among researchers asit does not fall victim to human bias as easily. Once the observations werebinned into certain grade groups, multinomial logistic regression was thenperformed on the data using the spectra as the predictors and the gradesas the response. This in theory provides the answer to the question ‘candegree of sickness be predicted using Raman spectroscopy?’ The resultingclassification tables can be found in Tables 4.3 and 4.4.As shown, the level of sickness can be predicted with a high level of accu-racy using the binning processes described above, yielding misclassificationrates of 2.15% and 3.15% respectively. This resulted in rates similar to ifthe original classes were kept. When using the original grades as groups, themultinomial logistic regression model still resulted in a misclassification rateof 2.94%.This section outlines how GBR-NMF can be effectively used as a methodfor removing unwanted data. It also displays that with the removal of spectra274.3. Raman Spectroscopic Cancer DataTable 4.3: Multinomial logistic regression predictions with respect to binningprocess number 1.PredictedTrue Low Medium HighLow: 0 - 0.1 2050 18 1Medium: 0.3 - 0.6 21 198 7High: 0.9 - 1 3 2 115Table 4.4: Multinomial logistic regression predictions with respect to binningprocess number 2.PredictedTrue Zero Low Medium HighZero: 0 1898 15 0 0Low: 0.1 - 0.3 31 264 1 0Medium: 0.5 - 0.6 6 18 62 1High: 0.9 - 1 2 0 2 116resembling media, level of sickness in mouse lung spectra can be predictedwith high levels of accuracy, regardless of how the grades are binned. Multiplebinning attempts were made and all of which provided similar results withrespect to classification accuracy using multinomial LASSO linear regression.Prediction of a single spectrum’s fibrosis grade can be interpreted in twoways; the probability of said cell being sick or the proportion of the sectionof the lung, which said cell belongs to, that is infected.28Chapter 5ConclusionA novel group and basis restricted nonnegative matrix factorization (GBR-NMF) algorithm was developed to input known bases and/or known groupswithin a nonnegative matrix factorization model. Updates were formulatedand then tested via simulations and applications to real data: specifically,digital facial images and Raman spectroscopy. GBR-NMF was used to specifyknown to exist a priori factors and groups in the models. Regarding the facialexpression data, GBR-NMF displayed significant improvements on standardNMF with respect to acquired factors from the added group constraint.Advancement over the standard NMF was also proven to be true concerningthe motivating data sets as the practicality of the algorithm greatly increase.Being able to specify known factors allows for spectra that may not havebeen found in the unconstrained form to be discovered.Future methodological work could consider adjustments made to theauxiliary matrix such as adjusting the off-diagonal elements to be of non-zerovalue. Further exploration of  in Equation 3.1 should also be considered.Namely, the development of new updating procedures in regards to differentnoise distributions and how the algorithm is affected when the assumption ofGaussian distributed error is violated. On the application side, the authorsand collaborators are currently investigating the unknown Raman spectrumacquired (fifth diagonal of Figure 4.5) from applying GBR-NMF to the datafrom Matthews et al. (2015) for potential biological importance in indicatingtumour radiosensitivity — such work has important ramifications for per-sonalized radiotherapy programs. Research regarding error distributions ofRaman spectroscopic data sets should also be investigated further, as theprocess of baseline subtraction and setting area under the curve of spectra toequal 1 could drastically affect this term.29BibliographyAbdi, H. and L. J. Williams (2010). Principal component analysis. Wileyinterdisciplinary reviews: computational statistics 2 (4), 433–459. 5Böhning, D. (1992). Multinomial logistic regression algorithm. Annals of theinstitute of Statistical Mathematics 44 (1), 197–200. 6Buciu, I. and I. Pitas (2004). Application of non-negative and local nonnegative matrix factorization to facial expression recognition. In Proceedingsof the 17th International Conference on Pattern Recognition, 2004. ICPR2004., Volume 1, pp. 288–291. IEEE. 18Butler, H. J., L. Ashton, B. Bird, G. Cinque, K. Curtis, J. Dorney,K. Esmonde-White, N. J. Fullwood, B. Gardner, P. L. Martin-Hirsch,et al. (2016). Using raman spectroscopy to characterize biological materials.Nature protocols 11 (4), 664. 2Dempster, A. P., N. M. Laird, and D. B. Rubin (1977). Maximum likelihoodfrom incomplete data via the em algorithm. Journal of the Royal StatisticalSociety: Series B (Methodological) 39 (1), 1–22. 13Deng, X., R. Ali-Adeeb, J. L. Andrews, P. Shreeves, J. J. Lum, A. Brolo, andA. Jirasek (2020). Monitor ionizing radiation-induced cellular responseswith raman spectroscopy, non-negative matrix factorization, and non-negative least squares. Applied Spectroscopy , 0003702820906221. 22Devarajan, K. (2008). Nonnegative matrix factorization: an analyticaland interpretive tool in computational biology. PLoS Comput Biol 4 (7),e1000029. 8Ding, C., X. He, and H. D. Simon (2005). On the equivalence of nonnegativematrix factorization and spectral clustering. In Proceedings of the 2005SIAM international conference on data mining, pp. 606–610. SIAM. 9, 11Feng, X., A. J. Moy, H. T. Nguyen, J. Zhang, M. C. Fox, K. R. Sebastian,J. S. Reichenberg, M. K. Markey, and J. W. Tunnell (2017). Raman activecomponents of skin cancer. Biomedical optics express 8 (6), 2835–2850. 230BibliographyFonti, V. and E. Belitser (2017). Feature selection using lasso. VU AmsterdamResearch Paper in Business Analytics 30, 1–25. 7Gray, R. M. (1990). Vector quantization. Readings in speech recognition 1 (2),75–100. 1Guillamet, D. and J. Vitria (2002). Classifying faces with nonnegative matrixfactorization. In Proc. 5th Catalan conference for artificial intelligence, pp.24–31. 18Harder, S. J., Q. Matthews, M. Isabelle, A. G. Brolo, J. J. Lum, and A. Jirasek(2015). A raman spectroscopic study of cell response to clinical doses ofionizing radiation. Applied spectroscopy 69 (2), 193–204. 21Harman, H. H. (1960). Modern factor analysis. 1Hartigan, J. A. and M. A. Wong (1979). Algorithm as 136: A k-meansclustering algorithm. Journal of the Royal Statistical Society. Series C(Applied Statistics) 28 (1), 100–108. 9Hastie, T., R. Tibshirani, and M. Wainwright (2015). Statistical learningwith sparsity: the lasso and generalizations. CRC press. 5Hilbe, J. M. (2009). Logistic regression models. CRC press. 7Hosmer Jr, D. W., S. Lemeshow, and R. X. Sturdivant (2013). Applied logisticregression, Volume 398. John Wiley & Sons. 7Huang, J., S. Ma, and C.-H. Zhang (2008). The iterated lasso for high-dimensional logistic regression. The University of Iowa, Department ofStatistics and Actuarial Sciences. 7Jolliffe, I. (2011). Principal component analysis. New York: Springer. 1, 5Jolliffe, I. T. (1986). Principal components in regression analysis. In Principalcomponent analysis, pp. 129–155. Springer. 5Jolliffe, I. T. (1990). Principal component analysis: a beginner’s guide—i.introduction and application. Weather 45 (10), 375–382. 4Jolliffe, I. T. and J. Cadima (2016). Principal component analysis: a reviewand recent developments. Philosophical Transactions of the Royal SocietyA: Mathematical, Physical and Engineering Sciences 374 (2065), 20150202.431BibliographyKleinbaum, D. G., K. Dietz, M. Gail, M. Klein, and M. Klein (2002). Logisticregression. New York: Springer. 7Kwon, S., S. Han, and S. Lee (2013). A small review and further studies onthe lasso. Journal of the Korean Data & Information Science Society 24,1077–1088. 5Lee, D. D. and H. S. Seung (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401 (6755), 788. 1, 7, 8, 10Lee, D. D. and H. S. Seung (2001). Algorithms for non-negative matrixfactorization. In Advances in neural information processing systems, pp.556–562. 8Li, T., C. Ding, and M. I. Jordan (2007). Solving consensus and semi-supervised clustering problems using nonnegative matrix factorization. InSeventh IEEE International Conference on Data Mining (ICDM 2007), pp.577–582. IEEE. 9Lin, C.-J. (2007). Projected gradient methods for nonnegative matrix factor-ization. Neural computation 19 (10), 2756–2779. 8Liu, H., Z. Wu, X. Li, D. Cai, and T. S. Huang (2011). Constrained nonneg-ative matrix factorization for image representation. IEEE Transactions onPattern Analysis and Machine Intelligence 34 (7), 1299–1311. 9, 10, 13, 14,17Matthews, Q., M. Isabelle, S. J. Harder, J. Smazynski, W. Beckham, A. G.Brolo, A. Jirasek, and J. J. Lum (2015). Radiation-induced glycogenaccumulation detected by single cell raman spectroscopy is associatedwith radioresistance that can be reversed by metformin. PloS one 10 (8),e0135356. 21, 22, 23, 29Matthews, Q., A. Jirasek, J. Lum, X. Duan, and A. G. Brolo (2010). Vari-ability in raman spectra of single human tumor cells cultured in vitro: cor-relation with cell cycle and culture confluency. Applied spectroscopy 64 (8),871–887. 21Paidi, S. K., P. M. Diaz, S. Dadgar, S. V. Jenkins, C. M. Quick, R. J.Griffin, R. P. Dings, N. Rajaram, and I. Barman (2019). Label-freeraman spectroscopy reveals signatures of radiation resistance in the tumormicroenvironment. Cancer research 79 (8), 2054–2064. 2132BibliographyPearson, K. (1901). Liii. on lines and planes of closest fit to systems of pointsin space. The London, Edinburgh, and Dublin Philosophical Magazine andJournal of Science 2 (11), 559–572. 4Pence, I. and A. Mahadevan-Jansen (2016). Clinical instrumentation andapplications of raman spectroscopy. Chemical Society Reviews 45 (7), 1958–1979. 2Psychological Image Collection at Sterling (2019). 2D face sets - pain expres-sions data at pics.stir.ac.uk. pics.stir.ac.uk. Accessed: 2019-09-12.18Shevade, S. K. and S. S. Keerthi (2003). A simple and efficient algorithmfor gene selection using sparse logistic regression. Bioinformatics 19 (17),2246–2253. 5Tibshirani, R. (1996). Regression shrinkage and selection via the lasso.Journal of the Royal Statistical Society: Series B (Methodological) 58 (1),267–288. 5Tibshirani, R. (1997). The lasso method for variable selection in the coxmodel. Statistics in medicine 16 (4), 385–395. 5Tibshirani, R. (2011). Regression shrinkage and selection via the lasso: aretrospective. Journal of the Royal Statistical Society: Series B (StatisticalMethodology) 73 (3), 273–282. 5Wagstaff, K., C. Cardie, S. Rogers, S. Schrödl, et al. (2001). Constrainedk-means clustering with background knowledge. In Icml, Volume 1, pp.577–584. 9Wang, Y.-X. and Y.-J. Zhang (2012). Nonnegative matrix factorization:A comprehensive review. IEEE Transactions on Knowledge and DataEngineering 25 (6), 1336–1353. 1, 8Wold, S., K. Esbensen, and P. Geladi (1987). Principal component analysis.Chemometrics and intelligent laboratory systems 2 (1-3), 37–52. 4Wu, C. J. (1983). On the convergence properties of the em algorithm. TheAnnals of statistics, 95–103. 13Yang, Z., Y. Zhang, Y. Xiang, W. Yan, and S. Xie (2018). Non-negativematrix factorization with dual constraints for image clustering. IEEETransactions on Systems, Man, and Cybernetics: Systems. 933BibliographyZdunek, R. and A. Cichocki (2007). Nonnegative matrix factorization withconstrained second-order optimization. Signal Processing 87 (8), 1904–1916.8Zhai, H. (2018). Variable selection via lasso with high-dimensional proteomicdata. Arts & Sciences Electronic Theses and Dissertations. 7Zhi, R., M. Flierl, Q. Ruan, and W. B. Kleijn (2010). Graph-preservingsparse nonnegative matrix factorization with application to facial expressionrecognition. IEEE Transactions on Systems, Man, and Cybernetics, PartB (Cybernetics) 41 (1), 38–52. 1834Appendix ADerivation of UpdatingAlgorithmDF =12|| X−WAS ||2∝ || X−WAS ||2= (√Tr((X−WAS)(X−WAS)T ))2= Tr((X−WAS)(X−WAS)T )= Tr((X−WAS)(XT − STATWT ))= Tr(XXT −XSTATWT −WASXT +WASSTATWT )= Tr(XXT )− Tr(XSTATWT )− Tr(WASXT )+Tr(WASSTATWT )= Tr(XXT )− Tr((XSTATWT )T )− Tr(WASXT )+Tr(WASSTATWT )= Tr(XXT )− Tr(WASXT )− Tr(WASXT )+Tr(WASSTATWT )= Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )As specified in section 3.2, this does not account for the nonnegativityconstraints. To do so, a Lagrange function L with multipliers αij , βij , andγij are required in order to constrain W, A, and S such that wij ≥ 0, aij ≥ 0and sij ≥ 0 respectively. This gives the following Lagrange function (withα = [αij ], β = [βij ], and γ = [γij ])L = Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )+Tr(αWT ) + Tr(βAT ) + Tr(γST )Now finding the derivatives of L with respect to W, A, and S we have35Appendix A. Derivation of Updating Algorithm∂L∂W=∂∂W(Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )+Tr(αWT ) + Tr(βAT ) + Tr(γST ))=∂∂WTr(XXT )− 2 ∂∂WTr(WASXT ) +∂∂WTr(WASSTATWT )+∂∂WTr(αWT ) +∂∂WTr(βAT ) +∂∂WTr(γST )= −2 ∂∂WTr(WASXT ) +∂∂WTr(WASSTATWT ) +∂∂WTr(αWT )= −2 ∂∂WTr(XSTATWT ) +∂∂WTr(WASSTATWT ) +∂∂WTr(αWT )= −2XSTAT + 2WASSTAT + α= 0∂L∂A=∂∂A(Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )+Tr(αWT ) + Tr(βAT ) + Tr(γST ))=∂∂ATr(XXT )− 2 ∂∂ATr(WASXT ) +∂∂ATr(WASSTATWT )+∂∂ATr(αWT ) +∂∂ATr(βAT ) +∂∂ATr(γST )= −2 ∂∂ATr(WASXT ) +∂∂ATr(WASSTATWT )∂∂ATr(βAT )= −2 ∂∂ATr(WTXSTAT ) +∂∂ATr(WTWASSTAT )∂∂ATr(βAT )= −2WTXST + 2WTWASST + β= 0∂L∂S=∂∂S(Tr(XXT )− 2Tr(WASXT ) + Tr(WASSTATWT )+Tr(αWT ) + Tr(βAT ) + Tr(γST ))=∂∂STr(XXT )− 2 ∂∂STr(WASXT ) +∂∂STr(WASSTATWT )+∂∂STr(αWT ) +∂∂STr(βAT ) +∂∂STr(γST )36Appendix A. Derivation of Updating Algorithm= −2 ∂∂STr(WASXT ) +∂∂STr(WASSTATWT ) +∂∂STr(γST )= −2 ∂∂STr(ATWTXST ) +∂∂STr(ATWTWASST ) +∂∂STr(γST )= −2ATWTX+ 2ATWTWAS+ γ= 0Each of these are set to equal zero as the goal is to find a local optimumof DF since it is unfeasible to find the global minimum. Applying the Karush-Kuhn-Tucker conditions, αijwij = 0, βijaij = 0, and γijsij = 0 to the threeequations results in the following(XSTAT )ijwij − (WASSTAT )ijwij = 0(WTXST )ijaij − (WTWASST )ijaij = 0 (A.1)(ATWTX)ijsij − (ATWTWAS)ijsij = 0which leads to the following updating ruleswij ← wij (XSTAT )ij(WASSTAT )ijaij ← aij (WTXST )ij(WTWASST )ijsij ← sij (ATWTX)ij(ATWTWAS)ij(A.2)37Appendix BProof of AlgorithmConvergenceAppendix B proves Lemmas 3.3 and 3.4 in order to fulfill the proofrequirements of Section 3.3.B.1 Proof of Lemma 3.3Proof. It is clear that G(a, a) = Faij (a). According to the definition of anauxiliary function, it only needs to be shown that G(a, atij) ≥ Faij (a). Thiscan be done using the Taylor series expansion of Faij :Faij (a) = Faij (atij) + F′aij (a− atij) +12F ′′aij (a− atij)2with F ′′aij being the second order derivative of Faij . We can show thatFaij = (∂D∂A)ij = (−2WTXST + 2WTWASST )ijF ′′aij = (∂2D∂A2)ij = 2(WTW)ii(SST )jjIt is now equivalent to prove(WTWASST )ijatij≥ 12F ′′ij = (WTW)ii(SST )jj (B.1)Now we have(WTWASST )ij =q∑`=1(WTWA)i`(SST )`j (B.2)≥ (WTWA)ij(SST )jj (B.3)≥q∑`=1(WTW)i`at`n(SST )nn (B.4)≥ atij(WTW)ii(SST )jj (B.5)38B.2. Proof of Lemma 3.4Which when rearranged proves Equation B.1.B.2 Proof of Lemma 3.4Proof. It is clear that G(s, s) = Fsij (s). According to the definition of anauxiliary function, it only needs to be shown that G(s, stij) ≥ Fsij (s). Thiscan be done using the Taylor series expansion of Fsij :Fsij (s) = Fsij (stij) + F′sij (s− stij) +12F ′′sij (s− stij)2with F ′′sij being the second order derivative of Fsij . We can show thatFsij = (∂D∂S)ij = (−2ATWTX+ 2ATWTWAS)ijF ′′sij = (∂2D∂S2)ij = 2(ATWTWAS)iiIt is now equivalent to prove(ATWTWAS)ijstij≥ 12F ′′ij = (ATWTWAS)ii (B.6)Now we have((ATWTWAS)ij =q∑`=1(ATWTWA)i`(S)`j (B.7)≥ (ATWTWAS)iistij (B.8)= stij(ATWTWAS)ii (B.9)which when rearranged proves Equation B.6.39
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Dimensionality reduction techniques with applications in Raman spectroscopy
Text