Multivariate Analysis of Fluorescence EEM Data

The dataset contains EEM measurements for OMP standards, natural waters, and wastewater samples collected in June and November. Each EEM is stored as a comma separated values (CSV) file with emission wavelengths as rows and excitation wavelengths as columns. Metadata files list sample identifiers, concentrations (conc_cip, conc_nap, conc_zol, in g/L), and descriptive labels such as natural water or wastewater influent and effluent. These metadata provide the supervised targets used in regression and classification.

The raw EEMs provided in the dataset were already corrected using standard fluorometric procedures from the source study, including instrument excitation and emission corrections, blank subtraction, and inner filter effect correction. Signals were normalized to the Raman peak area of ultrapure water, and scatter bands were removed and replaced with missing values.

Raw EEMs contain missing values around scatter bands and vary slightly in wavelength grids. All EEMs were standardized to a common grid by taking the intersection of wavelength axes across files. Missing values introduced by scatter removal showed up in the CSV files as NaN, which means “Not a Number”. These missing values were filled via smooth linear interpolation to provide a consistent feature layout for subsequent dimensionality reduction and model training.

Core preprocessing was performed using Pandas^[2], a tabular data library, with functions such as pandas.read_csv(), and NumPy^[3], a numerical computing library, with SciPy^[4] interpolation via scipy.interpolate.griddata().

Visualization of EEM Data

Heatmaps of raw and preprocessed EEMs were generated to visually inspect the effect of interpolation and scatter handling. The heatmaps of raw EEMs show scatter gaps and uneven intensity near the diagonal. After preprocessing, the filled regions are smooth and the overall fluorescence patterns remain intact. This comparison supports the choice of interpolation and indicates that the preprocessing step preserves the global structure required for downstream analysis, and the raw EEM heatmap with scatter gaps and the same sample after preprocessing are shown below. Plots are produced with Matplotlib^[5] using functions like matplotlib.pyplot.imshow().

Each preprocessed EEM is flattened into a feature vector and subjected to PCA. Principal component analysis is a linear technique that projects data into orthogonal (perpendicular) directions that capture the largest variance.

10 principal components were retained, which capture 99.93% of the total variance. This compresses each EEM into a small, stable feature set while preserving the dominant structure needed for prediction and visualization.

The PCA computation uses NumPy linear algebra with numpy.linalg.svd(). Here SVD stands for “Singular Value Decomposition” and this function is used because PCA is calculated as the SVD of the data matrix.

A permutation test is a non-parametric statistical method for testing whether two groups are significantly different regarding a chosen statistic. The core idea is to compare the observed statistic to a large sample of alternate statistics obtained by repeatedly shuffling group labels. Because the labels are randomly reassigned, this sample represents the null hypothesis that group membership does not matter. This approach avoids assumptions about normality or equal variances, which can be hard to justify for fluorescence data.

In this project, the statistic is the distance between the mean PCA vectors of the two groups. This statistic is simple and interpretable, and a large observed distance suggests that the groups occupy different regions of the reduced feature space, while a small distance suggests overlap. 5000 random relabelings were performed for each comparison and a p value was computed as the fraction of permuted distances greater than or equal to the observed distance. A small p value indicates that the observed separation is unlikely under random grouping. This test was applied to unspiked subsets of the data to isolate baseline source differences without spiking effects, which prevents concentration gradients from dominating the group comparison.

The test is implemented by concatenating samples from the two groups, shuffling the combined set, splitting it back into two groups of the original sizes, and recomputing the mean distance. This loop produces the null distribution of distances. The procedure is repeated for multiple group pairs, which enables the identification of which sources are clearly separated and which appear similar. A simple illustration of permutation testing is shown below.

Permutation testing is descriptive rather than predictive, and it complements the machine learning models. It provides a statistical check to confirm that the features contain real group differences before attempting to build models that exploit those differences to predict the metadata.

Because the test uses randomization, the p value can vary slightly between runs. A sufficiently large number of permutations was used to make this variability small relative to the effect sizes observed. The p values reported in the results section are therefore stable indicators of separation for the comparisons tested.

The permutation test also offers a way to interpret effect sizes. Larger observed mean distances indicate greater separation in the reduced space. When such distances coincide with small p values, the group differences are both statistically reliable and practically meaningful.

It is important to note that permutation testing does not explain why groups differ, only that they do. That is one reason to build predictive models, which help quantify how much information about spiking and source is contained in the EEM features.

A practical benefit of the permutation test is that it uses the same reduced features that feed the models. This means that the statistical test is directly tied to the features that the models rely on, so the two analyses reinforce each other rather than speaking about different spaces.

The permutation tests show that most source comparisons are significantly separated in PCA space, with the only non-significant comparison being between two natural water sources, likely due to similarity and small sample size. These results (shown below) support the interpretation that the reduced EEM features encode meaningful differences between sample origins.

The performance of multilayer perceptron and random forest models both depend on how certain parameters are chosen. These parameters that affect model performance, but must be determined before the model training process, are often called “hyperparameters” to distinguish them from any other parameters that are learned from data during training.

Grid search is a common way to find good hyperparameters. It works by evaluating a predefined set of parameter combinations and selects the best performing configuration according to a chosen metric. It provides a systematic way to explore tradeoffs such as model capacity, learning rate, or tree depth, and it makes the selection process transparent and repeatable. The alternative, manual tuning, is more subjective, less reproducible, and is likely to miss the best hyperparameter combinations.

Grid search has a computational cost because the number of combinations can grow quickly. To keep the search tractable, the focus was kept on ranges that are plausible given the dataset size and model complexity. For example, a moderate range of tree depths and hidden layer sizes were searched rather than extremely large values that would be unlikely to generalize.

An important consideration when training and evaluating models is which data the model should be trained on and which data it should be tested on. It is problematic to evaluate the performance of a model based on the data it has been trained on because models often memorize the training data and fail to generalize to new data. This is called “overfitting” and because of it, datasets are typically split into distinct training and testing sets where the model is only trained on the training set and tested on the testing set. This is also often extended with a “validation set” that is used as evaluation for iterative improvement of hyperparameters to preserve a testing set that has not been optimized for, even by hyperparameter selection. An alternative to this approach is to use 5-fold cross validation.

In 5-fold cross validation the dataset is split into five equal parts, called folds. The model is trained on four folds and tested on the remaining fold, and this is repeated so each fold serves as the test set once. Final performance is the average across the five test results. The concept of 5-fold cross validation is demonstrated in the video below.

This approach reduces dependence on any single train test split and provides a more reliable estimate of generalization performance. It is especially important when class sizes are imbalanced or when the dataset is modest in size, because a single split can overrepresent or underrepresent important groups.

We use stratified splitting when classification labels are involved so that each fold preserves the class distribution as much as possible. This keeps training and validation sets comparable and avoids misleading accuracy estimates. It also ensures that rare categories are present in both training and validation across folds. 5-fold cross validation is a balance between statistical reliability and computational cost. Using more folds would provide more training data in each split, but it would also increase computation time. For instance, for small datasets, sometimes “leave-one-out cross validation” is used where every datapoint is in its own fold. Five folds provides a stable estimate while keeping the workflow manageable.

A neural network is a learnable function that changes its shape to fit a set of input and output data. It is built from layers of weighted sums followed by non-linear activation functions, and the concept is shown here:

A multilayer perceptron (shown below) is the simplest form of neural network in which layers are fully connected and information flows forward from input to output.

The MLP is appropriate for fluorescence-based prediction because it can represent smooth non-linear relationships between principal component features and targets. It also supports multitask learning, which means a single model can predict concentrations and classifications simultaneously. This shared representation can improve efficiency and capture relationships between tasks, such as how spiking levels relate to source category.

The MLP is trained by minimizing a combined loss that includes regression error (MSELoss()) for concentration outputs and classification error (CrossEntropyLoss()) for categorical outputs. The model uses PyTorch^[7] components such as torch.nn.Linear() and torch.utils.data.DataLoader() to define layers and manage training batches.

The multitask setup requires balancing losses so that the regression and classification objectives both improve. The final architecture uses multiple hidden layers with non-linear activations to capture complex relationships while keeping the model small enough for the dataset size.

Grid search over the following MLP hyperparameters was performed: learning rate (2e-3, 1e-3, 7.5e-4, 5e-4, 1e-4), hidden layer sizes ("32,16", "64,32,16", "128,64,32", "256,128,64"), batch size (16, 32, 64), optimizer (stochastic gradient descent (SGD) or Adam), and whether or not we performed a log transformation of the concentration targets.

A log transform of concentration targets was considered to reduce the influence of large concentration values to help the model treat proportional errors more evenly across the concentration range and stabilize training when targets span orders of magnitude. While the best reported configuration did not rely on the log transform, the option remained useful for alternative settings. From an implementation perspective, the network uses small fully connected layers rather than convolutional layers because the input is already condensed into principal components (even if PCA had not been used, convolutional networks would not necessarily be suitable because the input EEM data is not translation invariant). This makes the model compact and efficient, which is important for repeated training during grid search and cross validation.

Conceptually, the MLP learns to reshape the input space so that patterns related to spiking and source become easier to separate. It does this by iteratively making small updates to its parameters in the direction that locally minimizes error on the training data. This process is called gradient descent, and when it is performed using only a subset of the training data for each update, it is called stochastic gradient descent. A visual illustration is shown below. This ability is central to why it outperforms the random forest on regression, where smooth relationships matter more (as seen in Section 4).

A random forest is an ensemble of decision trees that vote to make the final prediction for the random forest model. Each tree learns a set of simple “if-then” rules that partition the feature space, and the forest aggregates predictions across trees to reduce variance and improve generalization. For regression, the output is the average of tree predictions, and for classification, the output is the majority vote.

Decision Trees

A decision tree is a hierarchical model that recursively splits the feature space using simple rules, such as thresholds on a single feature. Each split aims to reduce prediction error for regression or increase class purity for classification. The resulting model is interpretable because each path from root to leaf corresponds to a clear sequence of rules. Decision trees can overfit when allowed to grow deep, but they capture non-linear decision boundaries without requiring feature scaling. A simple visual example of a decision tree is shown below.

The MLP was used because fluorescence patterns and concentration effects are expected to vary smoothly and non-linearly in the reduced space. A single network predicts three concentrations, a categorical source label, and a spiked or not spiked flag. This joint setup encourages shared representations across tasks. The trade-off is sensitivity to preprocessing and hyperparameters, and the risk of overfitting on small datasets.

The random forest was used because it is robust on tabular features and performs well with limited sample sizes. It also provides a strong baseline for classification. The downside is that it does not model smooth gradients as directly as the MLP, and separate models are needed for regression and classification.

Both model families are implemented in Python. The MLP uses the PyTorch deep learning library with functions such as torch.nn.Linear() and torch.utils.data.DataLoader(). The random forest models use scikit-learn, a machine learning library for Python, with functions such as sklearn.ensemble.RandomForestRegressor(), sklearn.ensemble.RandomForestClassifier(), and sklearn.model_selection.GridSearchCV().

Model selection was driven by the data and the research goals. The MLP is preferred when the main objective is accurate concentration regression, while the random forest provides a robust option for classification and a useful comparison baseline.

Because the two model families capture different kinds of relationships, their combined results strengthen the overall conclusions. When both models agree on classification outcomes, confidence that the source signal is robust is gained. When they differ on regression outcomes, insight is gained into the importance of smooth non-linear modeling for concentration prediction.

The choice to include both models also fits the educational goal of comparing different modeling philosophies. The MLP represents a flexible, data driven approach, while the random forest represents a structured ensemble approach with more built-in interpretability.