. Author manuscript; available in PMC: 2026 Sep 18.

Published in final edited form as: Comput Biol Med. 2025 Sep 18;197(Pt B):110959. doi: 10.1016/j.compbiomed.2025.110959

Tensor analysis of animal behavior by matricization and feature selection

Beichen Wang ^1,^*, Jiazhang Cai ^2,^¶, Luyang Fang ^2,^¶, Motokazu Tsujikawa ³, Ping Ma ², Yuk Fai Leung ^1,^4,^5,^6,^*

PMCID: PMC12510634 NIHMSID: NIHMS2109698 PMID: 40972230

Abstract

Contemporary neurobehavior research often collects multi-dimensional tensor (MDT) data containing time-series measurements for multiple features from multiple animals subjected to various perturbations. Proper analysis of the MDT data can reveal neural circuitries driving the behavior. However, many MDT analyses, such as tensor decomposition, may not yield results that are easy to interpret or directly compatible with standard multivariate analysis designed for 2-dimensional tensor (2DT) structures. To address this issue, the MDT data are transformed into 2DT by dimensionality reduction techniques, including matricization methods such as Index Construction and Feature Concatenation. Nonetheless, matricization may exclude key information or introduce spurious noise to multivariate analysis, so its impact on multivariate analysis remains elusive. Here, we demonstrated different matricization approaches and feature selection methods and evaluated their impacts on multivariate analysis performance using an MDT dataset of zebrafish visual-motor response collected from wildtypes and visually-impaired mutants. We matricized the MDT dataset using various Index Construction and Feature Concatenation methods, then identified informative 2DT features using the filter and embedded methods. To evaluate these feature-selection methods, we applied several classifiers to distinguish zebrafish of different genotypes and assessed their performances with cross-validation and holdout validation. We found that most classifiers performed the best using 2DT features matricized by Feature Concatenation and selected by the embedded method or union operation. The results also revealed unique behavioral differences between the wildtypes and mutants but not identified by multivariate analysis or MDT analysis. Our results demonstrate the utility of analyzing MDT behavioral data by matricization and feature selection.

Keywords: multi-dimensional tensor, matricization, animal behavior, feature engineer

1. Introduction

Animals display dynamic behaviors in response to environmental changes and stimuli, reflecting the underlying neural processing of their surroundings and experiences. Studying these animal behaviors facilitates the dissection of the neural processing at the cellular, circuitry, individual, and population levels [1-3]. The study approach can be observational [4] or experimental [5], wherein behavioral dynamics are recorded through one or more types of quantitative measurement [6,7]. Even when an experiment collects just one behavioral measurement in a time series, e.g., the displacement of animals, the resulting dataset contains multiple samples and time points. Such data are regarded as multivariate data (Fig 1A) and typically analyzed by multivariate analysis, which includes multivariate hypothesis testing (e.g., Hotelling’s T² Test [8], MANOVA [8], and PERMANOVA [9]) and machine learning (e.g., regression [10] and clustering [11,12]). As the variety of behavioral measurements increases with the advent of technology, a new issue emerges at the higher-level data structure, as illustrated in Fig 1B. In this data structure, even though each data point can be generally described by three factors or dimensions (color arrows in Fig 1B), each factor or dimension can contain multiple features (e.g., multiple distinct behavioral measurements, multiple samples with different genotypes and treatment conditions, or multiple time points). When this happens, this data structure can be generalized into a multi-dimensional tensor (MDT), also defined as multi-order tensors in tensor terminology [13-15] (see Section 2.1). In this study, we refer to MDT as tensors with more than two dimensions. Since the MDT structure contains more than two dimensions, it is incompatible with multivariate analysis. For example, MANOVA cannot directly use the MDT data as input. Hence, the MDT structure poses a challenge to data analysis in behavior studies.

Fig 1. — **(A)** The schematic structure of a time-series behavior data with one behavior measurement repeatedly collected from multiple samples (blue axis) over time (red axis). The resulting data structure is a matrix. Each brown rectangle represents a unique measurement per sample per time point. Its color intensity represents the magnitude of the measurement. **(B)** The schematic structure of a time-series behavior data with multiple behavior measurements (green axis) repeatedly collected from multiple samples (blue axis) over time (red axis). The resulting data structure is a three-dimensional array known as a three-dimensional tensor (3DT) or a third-order tensor, a type of multidimensional tensor (MDT). In this array, each “sheet” represents the data collected from one sample. The different color rectangles (brown, green, and blue) represent different behavior measurements per sample per time point, whereas their color intensities represent the magnitude of the measurement. **(C)** The schematic structure of MDT data when matricized by Feature Concatenation. The matricization process concatenates behavior measurements across a slice of the data cube in (B) at one time point, then iterates through all time points. As a result, the behavior measurement dimension is reduced, transforming the 3DT into a two-dimensional tensor or a matrix.

This challenge is partially resolved by transforming the MDT into a two-dimensional tensor (2DT) and reducing the number of dimensions. This transformation makes the transformed data structure compatible with multivariate analysis. The tensor dimensions can be reduced by many models, such as tensor component analysis based on CANDECOMP/PARAFAC (CP) tensor decomposition [13,16,17]. These models are generalizations of 2DT data dimensionality reduction. For instance, CP decomposition is conceptually similar to the principal component analysis (PCA) [17]. However, these models have not been widely applied to analyzing animal behaviors, possibly due to the difficulties in interpreting their mathematical representations, operations, and resulting outputs. The same interpretability issue occurs in PCA, where the resulting PCs are the mathematical representation of multiple original features. The abstract representation of the features makes it hard to interpret the significance and effect of the original features [18-20]. To mitigate this interpretability issue, many behavioral scientists reduce the tensor dimensions in their data with alternative approaches. One approach called matricization transforms features in MDT space into new features in a less complex 2DT space [13,21], which makes the resulting features more interpretable and compatible with multivariate analysis. Two matricization strategies are commonly used in behavior studies, which we refer to as 1) Index Construction and 2) Feature Concatenation. In Index Construction, a single index is constructed to represent all or a subset of the original features in the MDT [22]. For example, the average of a group of animal displacements can be calculated to reduce the sample dimension, assuming the samples in the group are homogenous [23-25]. In Feature Concatenation, all or a subset of features in a dimension of MDT are selected to form simpler 2DTs or matrices [26]. These matrices are then joined or concatenated, essentially flattening the complex data “cube” (Fig 1B) into a long data “sheet” (Fig 1C) [13,21,27] (See Section 2.1 for the technical details). For example, neural recordings across different trials can be joined and flattened onto other dimensions in MDT neural recording datasets [17,28]. Once the MDT data are reduced to 2DT data, they can be analyzed by multivariate analysis. However, both matricization strategies have drawbacks. For instance, Index Construction may exclude information from the original features that may be essential to describe the behavior dynamics, whereas Feature Concatenation may create too many transformed features and introduce additional noise to the downstream analysis. These issues may impact the performance of downstream multivariate analysis. Hence, both matricization strategies and their transformed features must be evaluated to select the best strategy and set of features that optimize the multivariate analysis performance.

This study aims to establish a feature selection framework for the matricization of MDT data in animal behavior studies and to compare the performance of different matricization and feature selection methods. We applied several matricization approaches to analyze an MDT dataset obtained from a zebrafish behavior assay. We empirically evaluated how different matricization and feature-selection strategies influence the performance of multivariate analysis, specifically in classification tasks. We found that the 2DT features selected by Feature Concatenation improved the classification performance the most. In Sections 2.1 and 2.2, we provide the relevant theoretical background on MDT and feature selection, respectively. Our study defined the optimal matricization strategy and feature selection methods for the example zebrafish behavior data set. It demonstrated a practical framework to matricize MDT behavior data and select important behavioral features to maximize the performance of downstream multivariate analysis.

2. Background

2.1. Multi-dimensional tensor

In this study, an MDT is defined as a tensor with more than two dimensions. For instance, a 3-dimensional tensor (3DT), $A$ , is represented in Fig 1B and can be denoted as $A ϵ R^{N \times M \times S}$ . In animal-behavior studies, $N$ represents the number of samples, $M$ represents the number of behavioral measurement types, and $S$ represents the number of time points. In other words, this 3DT can be considered a higher-order matrix or a matrix of matrices [13]. The 3DT can be transformed into a 2DT or a matrix using matricization methods such as Index Construction or Feature Concatenation. For example, to reduce the behavior-measurement dimension by Index Construction, a behavior index is constructed to represent all or a subset of $M$ features, reducing the number of features from $M$ to 1. The resulting 2DT, $A'$ , can be denoted as $A^{'} ϵ R^{N \times S}$ , equivalent to a matrix of $N$ rows by $S$ columns (as illustrated in Fig 1A). To reduce the behavior-measurement dimension by Feature Concatenation, the data is sliced by the number of features in a dimension to generate sub-matrices which are then joined to create one matrix. For example, in the case of Fig 1C, the concatenation essentially flattens the 3DT in Fig 1B to a 2DT, $A ″$ , denoted as $A^{″} ϵ R^{N \times M S}$ , equivalent to a matrix of $N$ rows by $M \times S$ columns. To reduce the time dimension, the time series can be transformed into symbolic representations [29,30] or behavior motifs [31,32]. In this case, the resulting matrix may contain $N$ rows by $M \times L$ columns, where $L$ denotes the reduced number of time points, i.e., $L < S$ .

2.2. Feature selection

After the matricization of an MDT, the transformed features in the 2DT space must be assessed to determine which ones effectively capture the behavioral dynamics. Only the informative features should be used in downstream multivariate analysis to optimize the analysis performance. The informative features can be selected by 1) filter method, 2) embedded method, and 3) wrapper method [33,34] (Table 1). The choice of method often relies on empirical experience and practical constraints [35].

Table 1. Three methods for feature selection.

Method	Approach	Examples	Advantages	Disadvantages
Filter	The filter method uses quantified metrics to rank and filter features without a machine-learning process [33,34,36].	P-value by statistical testing [37,38] Fold change [38] Correlation [39,40] Information Gain [33,41]	Low computational cost [33,34] Independent from the machine-learning processes [33,34]	Neglects dependencies among features [33,34]
Embedded	The embedded method relies on a built-in mechanism in a learning machine where learning and feature assessment are intertwined [33,34].	LASSO [42] Random Forest [43] Decision Tree [44]	Relatively low computation cost [33,34] Consider dependencies among features [33,34]	Feature importance depends on the mechanism and limitation of the classifier [33,34]
Wrapper	The wrapper method evaluates feature importance by its performance in a learning machine [33,34,36].	Sequential Forward Selection [45] Sequential Backward Elimination [45]	Consider dependencies among features [33,34] Superior performance for supervised machine learning [33,35,46]	High computation time for exhaustive search [33,34] Feature importance depends on the mechanism and limitation of the classifier [33,34]

Open in a new tab

To harness the power of all three methods and mitigate their drawbacks, some feature-selection workflows combine features selected by these methods through intersection or union operations, an approach known as the ensemble approach [47,48]. The ensemble approach via intersection has been shown to select features from complex medical datasets that achieve the best classification accuracy in empirical testing [48]. Therefore, it is possible that the ensemble approach can also be applied to complex behavior datasets. In this study, we demonstrated the feasibility of this approach by matricizing along the behavioral-measurement dimension. The wrapper method in the ensemble approach may be less compatible with behavior studies because the recordings may take hours or days, with a time resolution in seconds or milliseconds. There may be too many feature combinations to exhaustively iterate through under limited computational resources, especially for classifiers like neural networks, when applying the wrapper method. Therefore, we evaluated an ensemble framework that utilizes filter and embedded methods after reducing the dimensionality of an MDT. In our ensemble approach, the 2DT features were aggregated by intersection and union operations. The performance of both operations was evaluated by a classification task to analyze zebrafish visual behavior.

3. Methods

3.1. Zebrafish husbandry and breeding

The zebrafish used in this study were the AB strain obtained from the Zebrafish International Resource Center <https://zfin.org/ZDB-GENO-960809-7>. The visually-impaired transgenic line, Tg(rho:Hsa.RH1_Q344X), was generated and used in our previous studies and is referred to in this study as Q344X [23,49]. The animals were bred and maintained following standard procedures [50,51]. Prior to embryo collection, adults were sex-separated in breed tanks the night before breeding. They were mixed and bred from 9 a.m. to 11 a.m. the following day. All collected embryos were raised in the E3 medium at 28 °C in an incubator with the same light-dark cycle (14-hour light/ 10-hour dark) as in the fish facility. The medium was changed every day. During the process, unhealthy and dead embryos were discarded.

3.2. Visual-motor response assay

The MDT dataset used in this study was generated by the zebrafish visual-motor response (VMR) assay [52,53], a systematic approach to measure the visual startle (i.e., swimming behavior) of zebrafish larvae triggered by drastic light changes. In the VMR assay, the larvae are individually placed in a 96-well plate, and their visual startles are recorded by an infrared camera in a Zebrabox System (Viewpoint Life Sciences). Prior to the VMR assay, zebrafish larvae at 6 days post-fertilization were transferred to the E3 medium in Whatman UNIPlate square 96-well plates (VWR). The 96-well plates with larvae were placed in lightproof boxes for 14 hours of dark habituation inside the incubator. At 7 days post-fertilization, these 96-well plates were placed inside the Zebrabox for data collection. The larval VMR was triggered by two types of light stimulation: an onset (light-on) and an offset (light-off). The videos were recorded in the tracking mode, binning the behavioral measurements every second. For the light-on VMR, the protocol consisted of a 1-hour dark period and a 5-minute (300-second) light period. For the light-off VMR, the protocol consisted of a 1-hour dark period, a 1-hour light period, and a 5-minute (300-second) dark period. The VMR data used in this study were collected during our ongoing drug screening for retinitis pigmentosa using the Q344X mutant [23,49] to detect drug improvement in rod response [54]. The scotopic light intensity was fine-tuned to detect the maximal difference in the light-on and light-off VMR between wildtype and the Q344X mutant. The total irradiance for light stimulation over the visible spectrum was 0.00965 μW/cm² and 0.000389 μW/cm² for the light-on and light-off VMR, respectively. All experiments were conducted between 9 a.m. and 6 p.m. to minimize the effect of circadian rhythm on vision [23,52].

3.3. VMR data processing and analysis

3.3.1. General

All data analyses in this section were conducted with R statistical environment version 4.3.1 <https://www.r-project.org/>, except for data extraction from the VMR videos, which was conducted with ZebraLab version 3,22,3,77 (Viewpoint Life Sciences), and the multilayer perceptron classifier, which was conducted with Python version 3.9 <https://www.python.org/>.

3.3.2. VMR data structure

The larval activity in the recorded video was extracted and quantified into nine behavioral measurements (Fig 2A). First, the larval movement was categorized into three activity levels: inactivity ( $i n a$ ; speed < 0.6 cm/s), small activity ( $s m l$ ; speed ≥ 0.6 cm/s and < 1 cm/s), and large activity ( $l a r$ ; speed ≥ 1 cm/s). Speed was calculated based on every 5 frames of the video recording using the following equation, according to the specifications of Viewpoint Zebralab:

s p e e d = \frac{f \times \sum_{i = 1}^{4} d (t_{i + 1} - t_{i})}{4},

where $f$ denotes the frame rate of the camera; $t_{i}$ denotes the location of the animal center mass at Frame $i$ ; $d (t_{i + 1} - t_{i})$ denotes the distance of the animal center mass between Frame $t_{i + 1}$ and Frame $t_{i}$ . Second, within each activity category, three measurements were collected: duration of movement ( $d u r$ ), count of speed change ( $c t$ ), and distance of movement ( $d i s t$ ). Hence, nine behavioral measurements were taken per larva per second (Fig 2A). The complete definitions of these nine measurements are listed in the S1 Table.

Fig 2. — **(A)** The schematic of the 3DT of the VMR data. Zebrafish VMR was recorded using nine behavior measurements. Three types of behavior were recorded: duration ( $d u r$ ), count ( $c t$ ), and distance ( $d i s t$ ), each at three activity levels: inactivity ( $i n a; < 0.6$ cm/s), small activity ( $s m l; \geq 0.6$ cm/s and < 1 cm/s), and large activity ( $l a r; \geq 1$ cm/s). See Table S1 for the detailed definitions. Behavior measurements are represented by unique colors, whereas the activity levels are represented by color intensities, with the activity magnitude represented by bar sizes. These nine measurements were collected from each sample every second. The data will also contain multiple time points (s) and samples (n), resulting in 3DT data. **(B)** Index Construction by the Total Distance index. This index is the sum of the $d i s t$ measurements at three activity levels, i.e., $i n a$ , $s m l$ , and $l a r$ , per sample per second (left). This index is calculated for all samples (middle) over time (right). This results in a simpler data structure with one Total Distance index instead of nine behavior measurements. In the middle and right panels, each bar represents a Total Distance index for a sample in a second, whereas the grey tones represent the magnitude of the index. **(C)** Feature Concatenation by concatenating $d i s t$ measurements (Concat3Dist). The left panel represents the $i n a d i s t$ , $s m l d i s t$ , and $l a r d i s t$ from the original data (A). These features were joined/concatenated per sample per second (left). The concatenation is applied for all samples (middle) over time (right). This results in a simpler data structure arranged in a 2DT space. In the middle and right panels, each bar represents the three concatenated features in three grey tones for a sample in a second, whereas the grey tones represent the magnitude of the behavior measurement.

3.3.3. Matricizing the 3DT VMR data

Before matricization, the effect of confounding factors was removed from the VMR data to retain only the biological effect. Confounding factors included light intensity variations, batch effect, and baseline activity variations across different experiments [55]. These factors were removed using established methods (See S1 Appendix for details). Then, the normalized 3DT VMR data were matricized by Index Construction and Feature Concatenation. In Index Construction, an index was calculated to represent part or all of the behavioral measurements. For instance, a commonly used index in the VMR data is the Total Distance [23,56-58], the sum of distance ( $d i s t$ ) across all three activity levels (i.e., $i n a$ , $s m l$ , and $l a r$ ) for each sample (or larva) at each second for all time points (Fig 2B). In other words, the Total Distance is calculated as follows:

d i s t_{{n, s}} = i n a d i s t_{{n, s}} + s m l d i s t_{{n, s}} + l a r d i s t_{{n, s}},

where $d i s t_{{n, s}}$ denotes the Total Distance for sample $n$ at time point $s$ , with $n = 1, 2, \dots, N$ and $s = 1, 2, \dots, S$ . Here, $N$ denotes the total number of samples, $S$ denotes the total number of time points, and $i n a d i s t_{{n, s}}$ , $s m l d i s t_{{n, s}}$ , and $l a r d i s t_{{n, s}}$ denote the inactivity, small activity, and large activity distance, respectively, for sample $n$ at time point $s$ . Calculating the $d i s t$ behavior index reduces the number of behavioral measurements from 9 to 1 by merging 3 $d i s t$ measurements and discarding the remaining 6 measurements. This operation transforms the 3DT to a 2DT of $N \times S$ dimensions with $S$ features. In this study, we named this type of 2DT feature with $d i s t a n c e$ and the time point after the light change, separated by a dot (“.”). For example, $d i s t a n c e$ refers to the total-distance feature at 1 second after light change for each larva. We also tested two other indices, $L 1 Norm$ and $L 2 Norm$ , which utilize all 9 features instead of discarding some features that may contain unique information. In data science, $L 1$ and $L 2 Norms$ are used to reduce the dimensionality of the multivariate data [59]. Mathematically, $L 1$ and $L 2 Norms$ are defined as follows:

Given a vector:

x^{T} = [x_{1}, x_{2}, \dots, x_{n}],

where $x_{1}, x_{2}, \dots, x_{n}$ represent values for each dimension. The $L 1 Norm$ of $x$ is defined as:

L 1 N o r m = \sum_{i = 1}^{n} ∣ x_{i} ∣,

and the $L 2 Norm$ of $x$ is defined as:

L 2 N o r m = \sqrt{\sum_{i = 1}^{n} {∣ x_{i} ∣}^{2}} .

In our VMR data, the $n th$ larval movement at $s th$ second can be regarded as a vector of size $M$ . The corresponding $L 1$ and $L 2 Norms$ can be calculated as follows:

L 1 N o r m_{n s} = \sum_{i = 1}^{M} {∣ m_{i} ∣}_{n s}

L 2 N o r m_{n s} = \sqrt{\sum_{i = 1}^{M} {∣ m_{i} ∣}_{{n s}^{2}}}

where ${∣ m_{i} ∣}_{ns}$ denotes the $i th$ behavioral measurements of sample $n$ at time $s$ , with, $i = 1, 2, \dots M$ , $n = 1, 2, \dots, N$ , and $s = 1, 2, \dots, S$ . Since these behavioral measurements were in different units and scales, their values were first standardized to their corresponding means and standard deviations. After standardization, $L 1$ and $L 2 Norms$ were calculated as the 2DT features. The 2DT features were named following the same convention used in the Total Distance index. For example, $L 1 N o r m . 1$ refers to the $L 1 Norm$ feature at 1 second after light change for each larva.

In Feature Concatenation, a few or all behavioral measurements were concatenated per time point per sample to reduce the behavioral measurement dimension, as illustrated in Fig 1C. In this study, we reduced the behavior-measurement dimension by two Feature Concatenation approaches referred to as Concat3Dist and Concat9Meas. In Concat3Dist, we concatenated all three distance-related measurements (i.e., $i n a d i s t$ , $s m l d i s t$ , and $l a r d i s t$ ) for each larva at each second (Fig 2C). After joining each $d i s t$ measurement, the transformed matrix will have the dimension of $N \times (3 \times S)$ with $3 \times S$ features in the 2DT space, where $N$ denotes the size of the sample dimension and $S$ denotes the size of the time dimension. In Concat9Meas, we concatenated all nine behavioral measurements. After joining the measurements, the transformed matrix will have the dimension of $N \times (9 \times S)$ with $9 \times S$ features in the 2DT space. In both 2DT matricized by either Concat3Dist or Concat9Meas, each 2DT feature was named using the behavioral measurement and the time point after the light change, separated by a dot (“.”). For example, $i n a d i s t . 1$ refers to the feature representing the $i n a d i s t$ (inactivity distance) at 1 second after light change for each larva.

3.3.4. Feature selection

The 2DT features were selected by aggregating the results of the filter and embedded methods. For the filter method, two metrics were used: the p-value of the Wilcoxon rank sum test [60,61] and the average fold change ( $FC$ ) [62]. The p-value was calculated by running a Wilcoxon rank sum test between the wildtype and Q344X larvae using each feature and adjusted for multiple hypothesis testing using Holm’s method [63]. A feature was deemed important and selected if the adjusted p-value was less than 0.01. The $FC$ was calculated for each feature by the following formula:

F C = \frac{V_{i} (w i l d t y p e)}{V_{i} (Q 344 X)},

where $V_{i} (w i l d t y p e)$ and $V_{i} (Q 344 X)$ denote the average value of feature $i$ in all wildtype and Q344X larvae, respectively, with $i = 1, 2, \dots, S$ for features matricized by Index Construction and $i = 1, 2, \dots, M \times S$ for features matricized by Feature Concatenation. For both Index Construction and Feature Concatenation, $S$ denotes the number of time points; $M$ denotes the number of behavior-measurement types. A feature was selected if the absolute value of its $\log 2 (FC)$ was greater than 0.38 in light-on VMR (equivalent to a linear $FC > 1.3$ or < 0.77 [i.e., 1/1.3]) or 0.58 in light-off VMR (equivalent to a linear $FC > 1.5$ or < 0.67 [i.e., 1/1.5]). These cutoffs were empirically chosen to balance the stringency of selection and the number of selected features based on our previous observations that the wildtype larvae typically display a lower scotopic VMR under the light-on stimulus [54]. For the embedded method, Random Forest (RF) [43] was used. RF ranked the feature importance by the Mean Decrease of the Gini Index (MeanDecreaseGini) of each feature [64,65]. MeanDecreaseGini is calculated when RF randomly splits all features into random subsets and assesses information impurity for the features [64,66]. Finally, the features selected by the filter and embedded methods were aggregated by either intersection or union operation to create intersection and union sets.

3.3.5. Evaluation of feature selection and matricization performance

The performance of the different selected feature sets and matricization methods was evaluated by a classification task to distinguish the wildtype from the Q344X mutant. The following classifiers were used: Support Vector Machine (SVM), Naïve Bayes (NB), K-Nearest Neighbor (KNN), Decision Tree (DT), RF, Extreme Gradient Boosting (XGB), and Multilayer Perceptron (MLP). SVM, KNN, and NB are classic classifiers previously used for VMR analysis [67]. Tree-based methods, including DT, RF, and XGB, are considered superior in performance to the classic classifiers in other fields [68,69]. MLP is a type of neural network classifier [70] that is adaptive and capable of non-linear relationships in the data, achieving high accuracy and robustness [71]. Detailed information on these packages and their hyperparameter selection can be found in S2 Table. For each classifier, five feature sets were used: Full set, features selected by the filter method only (Filter), features selected by the embedded method only (Embedded), intersection of the features selected by the Filter and Embedded methods (Intersect), and union of the features selected by the Filter and Embedded methods (Union). The light-on and light-off VMR were analyzed separately.

The performance of the classifiers on each feature set was evaluated using an HV approach with 10-fold CV tuning (Fig 3). First, 80% of the dataset was allocated as the training set, and 20% of the dataset was allocated to the testing set. These training and testing sets were preprocessed separately to remove confounding factors prior to matricization. This training set was used in the 10-fold CV to tune each classifier’s hyperparameters and assess the performance of each feature set under different matricization methods. A range of hyperparameters per classifier was iterated using a grid search method [32,38], with the area under the receiver operating characteristic (AUROC) curve as the performance metric. The best hyperparameters for each classifier and feature set were selected based on the highest average AUROC score obtained in the CV. Next, the testing set (i.e., the remaining 20% of the data) was used in the HV to evaluate the generalizability of the feature sets and matricization methods. Each classifier was trained using the training set with the hyperparameters identified as optimal in the CV and then applied to the testing set. Classification performance on the testing set was measured using accuracy, sensitivity, specificity, precision, and the Kappa statistic [69,70].

4. Results

4.1. The 3DT VMR dataset

The 3DT behavior dataset used in this study consisted of a light-on and a light-off VMR dataset, each with 9 biological replications of the VMR collected from 24 wildtype and 24 Q344X larvae at 7 days post-fertilization. In the VMR protocol, the final session was a 5-minute (300-second) light-on or light-off period for triggering the light-on and light-off VMR, respectively. Hence, the analysis focused on this 5-minute period. The data structure followed the one outlined in Fig 2A and has 9 behavioral measurements ( $M = 9$ ), 432 samples ( $N = 432$ ), and 300 seconds after light change ( $S = 300$ ). Therefore, our input MDT has dimensions of 432×9×300. This data was matricized to 2DT data by Index Construction or Feature Concatenation. For instance, when a Total Distance index was constructed (Fig 2B), 3 distance measurements were reduced to 1 distance index, while the 6 other measurements were discarded. This operation transformed the VMR data from 3DT to 2DT format of 432×300 dimensions. The Total Distance is a commonly used index for this type of data, as it provides an intuitive way to visualize the behavioral and biological differences between sample types [23,56-58]. For example, the visually-impaired Q344X larvae lacked the visual startle immediately after light onset (Fig 4A) and light offset (Fig 4B), compared to the wildtype larvae. They also displayed a lower sustained response after light offset (Fig 4B). These intuitive differences allow for setting up simple rules and statistical tests to detect differences between genotypes, including using the fast response in the first few seconds to study visual startle [8,23,54]. In addition to Total Distance, the data can also be matricized by concatenating the 3 distance measurements along the time tensor, transforming the data to a matrix of 432×(3×300) or 432×900 in size (Fig 2C). This matrix can also be used for downstream multivariate analysis.

Fig 4. — **(A)** The light-on VMR of the wildtype (black trace) and Q344X (red trace) larvae. Each trace represents data collected from 9 biological replicates with 24 larvae per genotype per replicate. The plot shows the average Total Distance traveled per second by the two genotypes from 30 seconds before light onset at 0 to 300 seconds after light onset. The color ribbon on each trace represents the standard error of the mean (SEM). **(B)** The light-off VMR of the wildtype (black trace) and Q344X (red trace) larvae. Each trace represents data collected from 9 biological replicates with 24 larvae per genotype per replicate. The plot shows the average Total Distance traveled per second by the two genotypes from 30 seconds before light offset at 0 to 300 seconds after light offset. The color ribbon on each trace represents the SEM. The black and white boxes on each plot represent the light-off and light-on period.

4.2. Feature Interpretation

As illustrated in Section 4.1, the Total Distance index is easy to interpret and effectively describes the behavior dynamics of fish larvae. This ease of interpretability also applies to other matricization methods as well. For example, the other two indices in Index Construction, $L 1$ and $L 2$ , reflect the larval activity similarly to Total Distance—higher $L 1$ and $L 2$ values indicate greater larval activity, and vice versa. This is because $L 1$ and $L 2$ are scalars of the original behavioral measurements. In Feature Concatenation, Concat3Dist concatenates three behavioral measurements at each second: $i n a d i s t$ , $s m l d i s t$ , and $l a r d i s t$ . For instance, $l a r d i s t$ refers to the large activity distance, depicting the distance of larva movement at high speed (⩾ 1 cm/s) at each second. In other words, higher $l a r d i s t$ values indicate greater distances traveled at high speed. Concat9Meas concatenates all 9 behavioral measurements at each second: $i n d u r$ , $s m l d u r$ , $l a r d u r$ , $i n a c t$ , $s m l c t$ , $l a r c t$ , $i n a d i s t$ , $s m l d i s t$ , and $l a r d i s t$ . For instance, $l a r d u r$ refers to the large activity duration, depicting the duration of larval movement at high speed (⩾ 1 cm/s) at each second. In other words, higher $l a r d u r$ values indicate greater time spent traveling with fast and dashing maneuvers. Another measurement, $i n d u r$ , refers to the inactivity duration, depicting the duration of larvae movement with minimal speed (< 0.6 cm/s) at each second. In other words, higher $i n a d u r$ values indicate greater time spent traveling with slow and baseline maneuvers. Detailed definitions and indications of each measurement are provided in Table S1. Hence, by using matricization, the meaning of each behavioral measurement is retained in the Index Construction and Feature Concatenation features: when a larva has higher feature values, it indicates that the larva displayed a greater movement at the corresponding activity level at a given second.

In contrast, CP decomposition yields less interpretable features (S1 Fig). It projects data points from 3DT space to subspaces along each dimension. Although it revealed several unique patterns of rank scores across different ranks and factors, it is more challenging to interpret the meaning of higher or lower values of rank scores. For instance, in the light-on VMR data, the rank scores for Rank 1 of the Time factor showed an initial increase, followed by a subsequent decrease (S1 Fig A). Rank 1 of the Measurement factor showed that $i n a c t$ , count of speed change from inactivity (< 0.6 cm/s) to small activity (⩾ 0.6 cm/s and < 1 cm/s), showed the lowest score (S1 Fig B). Rank 1 of the Sample Factor showed rank score distributions were similar between wildtype and Q344X samples, except that wildtype samples showed slightly more dispersed rank scores than Q344X (S1 Fig C). However, these patterns and trends in rank scores do not explicitly show the absolute high or low behavior activities. In the Time factor, it is hard to interpret how much high- or low-speed maneuvers and how each larva contributed to the upward then downward trend; in the Measurement factor, it is hard to interpret at which time point and which larvae showed variations in $i n a c t$ ; in the Sample factor, it is hard to interpret at which time point or which measurements wildtype and Q344X showed variations in behavior. Hence, the result is not as directly interpretable as matricization. Overall, CP decomposition reveals unique patterns across ranks and factors, but it remains challenging to interpret the meaning of the rank scores and to understand how wildtype and Q344X differ in behavior across time points—an issue more effectively addressed by matricization.

4.3. Feature selection—filter method

The filter method was used to select 2DT features from the light-on and light-off VMR datasets matricized by Index Construction (Fig 5) and Feature Concatenation (Fig 6). Fig 5 shows the results of Index Construction by three indices: Total Distance, $L 1 Norm$ , and $L 2 Norm$ . In the light-on VMR dataset, the filter method selected only $d i s t a n c e . 1$ (i.e., the Total Distance in the first second) from the 300 features in the Total Distance index (Fig 5A). It did not select any features from the $L 1 Norm$ index (Fig 5B), but it selected $L 2 N o r m . 1$ from the 300 features in the $L 2 Norm$ index (Fig 5C). In the light-off VMR dataset, the filter method selected $d i s t a n c e . 2$ and $d i s t a n c e . 3$ from the 300 features in the Total Distance index. (Fig 5D). The filter method also selected $L 1 N o r m . 3$ from 300 features in the $L 1 Norm$ index (Fig 5E) and 3 of the 300 features in the $L 2 Norm$ index, namely, $L 2 N o r m . 2$ , $L 2 N o r m . 3$ , and $L 2 N o r m . 7$ (Fig 5F).

Fig 6 shows the feature-selection results of Feature Concatenation by Concat3Dist and Concat9Meas. In the light-on VMR data, the filter method did not select any feature from the 900 features (3 measurements x 300 s) in Concat3Dist (Fig 6A). It selected only $indur.1$ from the 2700 features (9 measurements x 300 s) in Concat9Meas (Fig 6B). In the light-off VMR, the filter method selected $lardist.2$ from the 900 features from Concat3Dist (Fig 6C). It also selected 37 of 2700 features in Concat9Meas (Fig 6D). Among these 548 features, $larct.2$ , $lardist.2$ , $larct.2$ and $inadur.7$ had the highest −log₁₀(p) and $∣ \log_{2} (FC) ∣$ values.

4.4. Feature selection—embedded method

Next, the embedded method was used to select 2DT features from the light-on and light-off VMR datasets matricized by Index Construction (Fig 7) and Feature Concatenation (Fig 8). Fig 7 shows the results of Index Construction by three indices: Total Distance, $L 1 Norm$ , and $L 2 Norm$ . In the light-on VMR dataset, the embedded method selected $d i s t a n c e . 1$ from the 300 features in the Total Distance Index (Fig 7A). The method also selected $L 1 N o r m . 1$ and $L 1 N o r m . 41$ from the 300 features in the $L 1 Norm$ index (Fig 7B). It also selected 5 features from the 300 features in the $L 2 Norm$ index (Fig 7C), namely, $L 2 N o r m . 1$ , $L 2 N o r m . 41$ , $L 2 N o r m . 69$ , $L 2 N o r m . 70$ , and $L 2 N o r m . 96$ . In the light-off VMR dataset, the embedded method selected 5 from the 300 features from the Total Distance index (Fig 7D), namely, $d i s t a n c e . 2$ , $d i s t a n c e . 3$ , $d i s t a n c e . 6$ , $d i s t a n c e . 7$ , and $d i s t a n c e . 53$ . The embedded method also selected 9 from the 300 features in the $L 1 Norm$ index (Fig 7E). Among these 9 features, $L 1 N o r m . 2$ , $L 1 N o r m . 3$ , and $L 1 N o r m . 7$ had the highest values of MeanDecreaseGini. In addition, the embedded method selected 9 from the 300 features in the $L 2 Norm$ index (Fig 7F). Among these 9 features, $L 2 N o r m . 2$ , $L 2 N o r m . 3$ , and $L 2 N o r m . 7$ had the highest MeanDecreaseGini values.

Fig 8 shows the feature-selection results of Feature Concatenation by two approaches: Concat3Dist and Concat9Meas. In the light-on VMR data, the embedded method selected 3 from the 900 features in Concat3Dist (Fig 8A), namely, $i n a d i s t . 1$ , $s m l d i s t . 1$ , and $i n a d i s t . 131$ . The embedded method also selected 5 from the 2700 features in Concat9Meas (Fig 8B), namely, $s m l c t . 1$ , $s m l d u r . 1$ , $i n a d i s t . 1$ , $i n a d u r . 1$ , and $i n a c t . 1$ . In the light-off VMR data, the embedded method selected 12 from the 900 features in Concat3Dist (Fig 8C). Among these 12 features, $l a r d i s t . 2$ , $s m l d i s t . 2$ and $l a r d i s t . 3$ had the highest MeanDecreaseGini values. The embedded method also selected 13 from 2700 features in Concat9Meas (Fig 8D). Among these 13 features, $l a r d i s t . 2$ , $l a r c t . 2$ , and $l a r d u r . 2$ had the highest MeanDecreaseGini values.

4.5. Aggregation of selected 2DT features and performance evaluation

Table 2 lists the number of 2DT features selected by the filter and embedded methods, along with their intersection and union sets. It also includes the feature reduction rates for each selection method. In general, the number of features was reduced most substantially when the features were selected by the intersection operation, regardless of the dataset (i.e., light-on or light-off). Even though this operation had the highest reduction rate (100%) for the corresponding matricization, it would not provide any feature for downstream analyses. Hence, the operation with the minimum non-zero number of features was used in the next analysis.

Table 2. The number of features selected by different approaches.

Feature		Index Construction						Feature Concatenation
		Total Distance (3 measurements)		$L 1 Norm$ (9 measurements)		$L 2 Norm$ (9 measurements)		Concat3Dist (3 measurements)		Concat9Meas (9 measurements)
		Number	Reduction Rate (%)	Number	Reduction Rate (%)	Number	Reduction Rate (%)	Number	Reduction Rate (%)	Number	Reduction Rate (%)
Light-on	Full	300	0.00	300	0	300	0	1800	0	2700	0.00
	Filter	1	99.67	0	100	1	99.67	0	100	1	99.96
	Embedded	1	99.67	2	99.33	5	98.33	3	99.83	5	99.81
	Intersection	1	99.67	0	100	1	99.67	0	100	1	99.96
	Union	1	99.67	2	99.33	5	98.33	3	99.83	5	99.81
Light-off	Full	300	0.00	300	0.00	300	0.00	1800	0.00	2700	0.00
	Filter	2	99.33	1	99.67	3	99.00	1	99.94	37	98.63
	Embedded	5	98.33	9	97.00	9	97.00	12	99.33	13	99.52
	Intersection	2	99.33	1	99.67	3	99.00	1	99.94	10	99.63
	Union	5	98.33	9	97.00	9	97.00	12	99.33	40	98.52

Open in a new tab

The Reduction Rate reflects the reduction of the feature numbers after matricization compared to the full feature set, expressed in percentage. For each light condition and matricization method, the bold numbers highlight the biggest reduction (excluding 100%) across all feature sets.

To evaluate the impact of feature selection on multivariate analysis performance, the selected 2DT features were used in a classification task to distinguish the wildtype and Q344X larvae by SVM, NB, KNN, DT, RF, XGB, and MLP. A CV was conducted on the training set with 80% of the entire data, using the features selected from three constructed indices (Total Distance, $L 1 Norm$ , and $L 2 Norm$ ) and two concatenated features (Concat3Dist and Concat9Meas). The resulting average AUROC for each classifier is listed in Table 3, and the optimal hyperparameters are shown in S3 Table. The performance was then evaluated by comparing the highest AUROC achieved by each classifier using features from data matricized by different approaches.

Table 3.

Average AUROC in CV using feature sets from Index Construction and Feature Concatenation.

	Classifier	Feature Set	Index Construction			Feature Concatenation
	Classifier	Feature Set	Total Distance (3 measurements)	$L 1 Norm$ (9 measurements)	$L 2 Norm$ (9 measurements)	Concat3Dist (3 measurements)	Concat9Meas (9 measurem ents)
L	SVM	Full	0.5928	0.6020	0.6088	0.6002	0.6367
		Filter	0.6647^*	NA	0.6200	NA	0.6375
		Embedded	0.6647^*	0.6087	0.6646	0.6731^*	0.6244
		Intersection	0.6647^*	NA	0.6200	NA	0.6375
		Union	0.6647^*	0.6087	0.6646	0.6731^*	0.6244
	NB	Full	0.5695	0.5693	0.5171	0.6038	0.5749
		Filter	0.6730^*	NA	0.6266	NA	0.6397
		Embedded	0.6730^*	0.6023	0.6057	0.6158	0.6541^*
		Intersection	0.6730^*	NA	0.6266	NA	0.6397
		Union	0.6730^*	0.6023	0.6057	0.6158	0.6541^*
	KNN	Full	0.5769	0.5624	0.5413	0.5807	0.5553
		Filter	0.6612	NA	0.6235	NA	0.6307
		Embedded	0.6637	0.6246	0.6818^*	0.6976	0.7329^*
		Intersection	0.6672	NA	0.6215	NA	0.6333
		Union	0.6612	0.6246	0.6818^*	0.6976	0.7329^*
	DT	Full	0.6574^*	0.6170	0.6108	0.6390	0.6438
		Filter	0.6555	NA	0.5976	NA	0.6407
		Embedded	0.6555	0.6443	0.6361	0.6630^*	0.6340
		Intersection	0.6555	NA	0.5976	NA	0.6407
		Union	0.6555	0.6443	0.6361	0.6630^*	0.6340
	RF	Full	0.6726	0.6536	0.6775	0.7469	0.8063^*
		Filter	0.5857	NA	0.6211	NA	0.5777
		Embedded	0.5857	0.6446	0.6912	0.6764	0.7067
		Intersection	0.5857	NA	0.6211	NA	0.5777
		Union	0.5857	0.6514	0.6996^*	0.6679	0.7063
	XGB	Full	0.6644	0.6234	0.6287	0.6935	0.6960
		Filter	0.6494	NA	0.6258	NA	0.6378
		Embedded	0.6494	0.6451	0.6968^*	0.7069^*	0.7069^*
		Intersection	0.6494	NA	0.6258	NA	0.6378
		Union	0.6494	0.6451	0.6968^*	0.7069^*	0.6960
	MLP	Full	0.5310	0.5549	0.5471	0.5652	0.6110
		Filter	0.6700	NA	0.6366	NA	0.6408
		Embedded	0.6693	0.6179	0.6509	0.6894	0.7184^*
		Intersection	0.6643	NA	0.6301	NA	0.6421
		Union	0.6760^*	0.6175	0.6470	0.6768	0.7131
Light-off	SVM	Full	0.6926	0.5933	0.5850	0.6313	0.6672
		Filter	0.7386	0.6577	0.7167	0.7043	0.7076
		Embedded	0.7414^*	0.7307	0.7130	0.7511^*	0.7493
		Intersection	0.7386	0.6577	0.7167	0.7043	0.7277
		Union	0.7414^*	0.7307	0.7130	0.7511^*	0.7110
	NB	Full	0.6389	0.6107	0.6132	0.6558	0.6403
		Filter	0.7352	0.6518	0.6988	0.6863	0.7407
		Embedded	0.7442^*	0.7277	0.6884	0.7572^*	0.7155
		Intersection	0.7352	0.6518	0.6988	0.6863	0.7092
		Union	0.7442^*	0.7277	0.6884	0.7572^*	0.7387
	KNN	Full	0.5726	0.5651	0.5601	0.4984	0.5216
		Filter	0.7563^*	0.6556	0.7340	0.7051	0.7368
		Embedded	0.7463	0.7308	0.7230	0.7791^*	0.7392
		Intersection	0.7395	0.6591	0.7263	0.7050	0.7268
		Union	0.7451	0.7277	0.7230	0.7781	0.7398
	DT	Full	0.7360	0.6872	0.6932	0.7276	0.6942
		Filter	0.7336	0.6768	0.7456	0.7263	0.7742
		Embedded	0.7596^*	0.7497	0.7470	0.7996^*	0.6564
		Intersection	0.7336	0.6768	0.7456	0.7263	0.6581
		Union	0.7596^*	0.7497	0.7470	0.7996^*	0.7907
	RF	Full	0.7268	0.7361	0.7222	0.7334	0.7633
		Filter	0.6808	0.6556	0.7102	0.6585	0.7505
		Embedded	0.7047	0.7382	0.7573^*	0.7637^*	0.7308
		Intersection	0.6843	0.6556	0.7028	0.6585	0.7223
		Union	0.7028	0.7375	0.7561	0.7534	0.7479
	XGB	Full	0.6988	0.7084	0.7100	0.6887	0.6946
		Filter	0.7114	0.6801	0.6877	0.6950	0.7428
		Embedded	0.7521	0.7380	0.7559^*	0.7443	0.7301
		Intersection	0.7114	0.6801	0.6877	0.6950	0.7226
		Union	0.7521	0.7473	0.7559^*	0.7486^*	0.7293
	MLP	Full	0.6708	0.6467	0.6267	0.6901	0.6769
		Filter	0.7379	0.6621	0.7360	0.6997	0.7668
		Embedded	0.7577^*	0.7086	0.7039	0.7783	0.7342
		Intersection	0.7341	0.6758	0.7348	0.6997	0.7206
		Union	0.7544	0.7049	0.7205	0.7924^*	0.7126

Open in a new tab

The numbers in boldface indicate the highest average AUROC across a row. In other words, it is the highest value for a specific classifier using a specific feature set across all matricization methods. The numbers highlighted by an asterisk indicate the highest average AUROC across different feature sets per matricization strategy (Index Construction or Feature Concatenation) per classifier. NAs indicate that condition does not have any selected feature for the calculation.

First, for the 2DT features selected from the light-on dataset matricized by Index Construction, 3 of the 7 classifiers achieved the highest average AUROC using the features selected from the Total Distance index (Table 3), namely, SVM (0.6647), NB (0.6730), and MLP (0.6760); 3 of the 7 classifiers achieved the highest average AUROC using the features selected from the $L 2 Norm$ index, namely, KNN (0.6818), RF (0.6996), and XGB (0.6968); 1 of the 7 classifiers achieved the highest average AUROC using the full feature set from the Total Distance, namely, DT (0.6574). Among the 6 classifiers that achieved the highest AUROC using selected features, 2 classifiers used either the embedded method or the union operation, namely, KNN (0.6818) and XGB (0.6968); RF (0.6996) and MLP (0.6760) used the union operation; 2 classifiers used any selected feature set to achieve the same highest AUROC, namely, SVM (0.6647) and NB (0.6730). For the 2DT features selected from the light-off dataset matricized by Index Construction, 5 of the 7 classifiers achieved the highest AUROC using the features selected from the Total Distance index. Among these 5 classifiers, SVM (0.7414), NB (0.7442), and DT (0.7596) used the features selected by the embedded method or the union operation; KNN (0.7563) used the features selected by the filter method; MLP (0.7577) used the features selected by the embedded method. Two classifiers achieved the highest AUROC using the features selected from the $L 2 Norm$ index: RF (0.7573) used the features selected by the embedded method, and XGB (0.7559) used the features selected by the embedded method or the union operation. Overall, most classifiers achieved the highest AUROC when using features selected from the data matricized by the Total Distance index.

Second, for the 2DT features selected from the light-on dataset matricized by Feature Concatenation, 3 of the 7 classifiers achieved the highest AUROC using the features selected from Concat9Meas rather than from Concat3Dist, namely NB (0.6541), KNN (0.7329), and MLP (0.7184); 2 of the 7 classifiers achieved the highest AUROC using the features selected from Concat3Dist rather than from Concat9Meas, namely SVM (0.6731) and DT (0.6630); 1 classifier, XGB (0.7069), achieved the highest AUROC using selected feature from either Concat3Dist and Concat9Meas; 1 classifier, RF (0.8063), achieved the highest AUROC using the full feature set from Concat9Meas. Among the 3 classifiers that achieved the highest AUROC using features selected from Concat9Meas, NB (0.6541) and KNN (0.7329) used features selected by the embedded method or the union operation, while MLP (0.7184) used the features selected by the embedded method. Two classifiers achieved the highest AUROC using features selected from Concat3Dist: SVM (0.6731) and DT (0.6630), both of which used the features selected by the embedded method or the union operation. One classifier, XGB (0.7069), used features selected by the embedded method or the union operation in Concat3Dist, and features selected by the embedded method in Concat9Meas. For the 2DT features selected from the light-off dataset matricized by Feature Concatenation, all 7 classifiers achieved higher AUROC using the features selected from Concat3Dist than those selected from Concat9Meas Among them, SVM (0.7511), NB (0.7572), and DT (0.7996) used the features selected by the embedded method or union operation; KNN (0.7791) and RF (0.7637) used the features selected by the embedded method; XGB (0.7486) and MLP (0.7924) used the features selected by the union operation. Overall, in the light-on dataset, most classifiers achieved the highest AUROC when using features selected from the data matricized by Concat9Meas; in the light-off dataset, most classifiers achieved the highest AUROC when using features selected from the data matricized by Concat3Dist.

Finally, to evaluate the performance across both matricization methods, we compared the AUROC of different classifiers using different feature sets selected from the data matricized either by Index Construction or Feature Concatenation. In the light-on dataset, 6 of the 7 classifiers had higher AUROCs using the features selected from Feature Concatenation than those from Index Construction, regardless of the feature set. Among the 6 classifiers, KNN (0.7329), RF (0.8063), and MLP (0.7184) used features selected from Concat9Meas; SVM (0.6731) and DT (0.6630) used features selected from Concat3Dist; XGB (0.7069) used features selected from either Concat9Meas or Concat3Dist. The only exception was NB (0.6730), which achieved a higher AUROC using features selected from the TotalDistance index. In the light-off dataset, 6 out of 7 classifiers had higher AUROCs using the features selected from Concat3Dist. The only exception was XGB (0.7559), which used features selected from the $L 2 Norm$ index. Overall, most classifiers achieved higher average AUROCs when using features selected from the data matricized by Feature Concatenation rather than Index Construction.

4.6. Generalization of the matricization approach and selected feature set to new data

After identifying the best-performing matricization approach and the resulting selected 2DT feature set in the CV, we then evaluated the extent to which their performances generalized to new data by HV (Fig 3). We first trained all classifiers using the training dataset (80% of the entire data) and all feature sets from all matricization methods, applied the optimal hyperparameters found in CV (S3 Table), and assessed their performances on the testing dataset (20% of the entire data) (S4 Table).

For Index Construction in the light-on dataset, 3 of the 7 classifiers achieved the highest AUROC using features selected from the Total Distance index, including KNN (0.6705), DT (0.5568), and XGB (0.6477); one classifier, SVM (0.6136), achieved the highest AUROC using features selected from the $L 1 Norm$ index; one classifier, NB (0.7273), achieved the highest AUROC using the full feature set from the $L 1 Norm$ index; one classifier, MLP (0.6705), achieved the highest AUROC using the full feature set from the $L 2 Norm$ index; and one classifier, RF (0.6023), achieved the highest AUROC using the full feature set from either the $L 1$ or $L 2 Norm$ indices. In the light-off dataset, 4 of the 7 classifiers achieved the highest AUROC using features selected from the Total Distance index, including SVM (0.7614), KNN (0.7727), DT (0.7273), and MLP (0.7727); one classifier, NB (0.6818), achieved the highest AUROC using features selected from the $L 1 Norm$ index; one classifier, XGB (0.7159), achieved the highest AUROC using features selected from the $L 2 Norm$ index; and one classifier, RF (0.7614), achieved the highest AUROC using the full feature set from the $L 1 Norm$ index. Overall, in HV, the classifiers varied the index that achieved the highest AUROC. Nonetheless, the greatest number of classifiers achieved the highest AUROC using features from the Total Distance index. These findings in HV differ from those in CV, where most classifiers achieved the highest AUROC using the Total Distance index for Index Construction.

For Feature Concatenation in the light-on dataset, 3 of the 7 classifiers achieved the highest AUROC using features selected from Concat3Dist, including KNN (0.6818), DT (0.6364), and MLP (0.6364); 3 of the 7 classifiers achieved the highest AUROC using features selected from Concat9Meas, including SVM (0.6136), NB (0.7368), and XGB (0.6818); one classifier, RF (0.7045), achieved the highest AUROC using the full feature set from Concat3Dist. In the light-off dataset, 6 of the 7 classifiers achieved the highest AUROC using features selected from Concat3Dist, including SVM (0.7841), NB (0.7273), KNN (0.8182), DT (0.7727), XGB (0.7273), and MLP (0.7841); One classifier, RF (0.7273), achieved the highest AUROC using features selected from Concat9Meas. Overall, these findings in HV are consistent with those in CV: the classifiers achieved the highest AUROC using features selected from Concat3Dist or Concat9Meas in the light-on dataset in both CV and HV; classifiers achieved the highest AUROC using features selected from Concat3Dist in the light-off dataset in both HV and CV.

We then identified the best-performing matricization strategy—Index Construction or Feature Concatenation—by comparing the highest AUROC of HV achieved by each classifier under each approach (Table 4). In the light-on dataset, 4 of the 7 classifiers achieved the highest AUROC using features selected from the dataset matricized by Feature Concatenation. Among these 4 classifiers, NB (0.7386), KNN (0.6818), and DT (0.6364) used features selected by the embedded method or the union operation; XGB (0.6818) used features selected by the union operation. One classifier, SVM (0.6136), achieved the highest AUROCs using features selected by either Index Construction or Feature Concatenation. One classifier, RF (0.7045), achieved the highest AUROC using the full feature set from Feature Concatenation. One classifier, MLP (0.6705), achieved the highest AUROC using the full feature set from Index Construction. In the light-off testing dataset, 6 of the 7 classifiers achieved the highest AUROC using the features selected from the data matricized by Feature Concatenation. Among these 6 classifiers, SVM (0.7841), NB (0.7273), DT (0.7727), XGB (0.7273), and MLP (0.7841) used the features selected by the filter method or the intersection operation; KNN (0.8182) used the features selected by the union operation. One classifier, RF, achieved the highest AUROC using the full feature set from Index Construction. Overall, most classifiers achieved the highest AUROC using the features selected by Feature Concatenation.

Table 4. The highest AUROC of HV in Index Construction and Feature Concatenation.

	Classifiers	Index Construction					Feature Concatenation
	Classifiers	Full	Filter	Embedded	Intersection	Union	Full	Filter	Embedded	Intersection	Union
Light-on	SVM	0.4318	NA	0.6136	NA	0.6136	0.5000	0.6136	0.5909	0.6136	0.5909
	NB	0.7273	NA	0.4432	NA	0.4432	0.5227	0.5795	0.7386	0.5795	0.7386
	KNN	0.5341	0.6591	0.6705	0.6250	0.6591	0.4659	NA	0.6818	NA	0.6818
	DT	0.5455	0.5568	0.5568	0.5568	0.5568	0.5795	NA	0.6364	NA	0.6364
	RF	0.6023	NA	0.5795	NA	0.5682	0.7045	NA	0.6477	NA	0.6136
	XGB	0.5795	0.6477	0.6477	0.6477	0.6477	0.6591	0.4659	0.5227	0.4659	0.6818
	MLP	0.6705	0.4432	0.5114	0.5000	0.5114	0.4091	NA	0.6364	NA	0.6364
Light-off	SVM	0.5682	0.7614	0.7273	0.7614	0.7273	0.5682	0.7841	0.7614	0.7841	0.7614
	NB	0.6591	0.6364	0.6818	0.6364	0.6818	0.5227	0.7273	0.5909	0.7273	0.5909
	KNN	0.6023	0.7727	0.7500	0.7614	0.7386	0.5000	0.7614	0.8068	0.7500	0.8182
	DT	0.6136	0.6932	0.7273	0.6932	0.7273	0.6136	0.7727	0.6591	0.7727	0.6136
	RF	0.7614	0.6023	0.7273	0.6250	0.7159	0.5795	0.7045	0.6364	0.7045	0.6477
	XGB	0.5568	0.6932	0.7159	0.6932	0.7159	0.5455	0.7273	0.6932	0.7273	0.7159
	MLP	0.7159	0.7727	0.7273	0.7727	0.7273	0.5227	0.7841	0.7386	0.7841	0.7273

Open in a new tab

The values in the table are the AUROCs achieved by each classifier on the testing dataset. The AUROCs shown here are the highest AUROC per matricization method for each classifier. The full performance results are provided in S4 Table. The numbers highlighted in bold typeface indicate the highest value in each row, i.e. the highest AUROC achieve by each classifier across all feature sets and matricization methods. NAs indicate that the corresponding condition did not yield a selected feature for the calculation. In light-on dataset, for Index Construction, KNN, DT and RF used features from the Total Distance index; SVM, NB, and XGB used features from the $L 1 Norm$ index; MLP used features from the $L 2$ index. For Feature Concatenation, SVM, NB, and RF used features from Concat9Meas; KNN, DT, XGB, and MLP used features from Concat3Dist. In the light-off dataset, for the Index Construction, SVM, KNN, DT, and MLP used the Total Distance index; NB and XGB used the $L 1 Norm$ Nom index; RF used the $L 2 Norm$ index. For the feature Concatenation, all except XGB used Concat3Dist; XGB used Concat9Meas

4.7. Performance evaluation of the matricization approach and selected features by other performance metrics

In addition to AUROC, other metrics can be used to evaluate performance. These metrics were used to evaluate the extent to which the performance conclusions drawn from the AUROC metric could be generalized to other performance metrics. Table 5 lists the values of various performance metrics for different classifiers. These values were calculated using the same best-performing feature set that achieved the highest AUROC in HV for the corresponding classifier shown in Table 4.

Table 5.

Performance by different metrics in HV using the best-performing feature set.

	Classifiers	AUROC	Accuracy	Sensitivity	Specificity	Precision	Kappa
Light-On	SVM	0.6136	0.6136	0.5227	0.7045	0.6389	0.2273
	NB	0.7386	0.7386	0.7273	0.7500	0.7442	0.4773
	KNN	0.6818	0.6818	0.7045	0.6591	0.6739	0.3636
	DT	0.6136	0.6136	0.6591	0.5682	0.6042	0.2273
	RF	0.6932	0.6932	0.9545	0.4318	0.6269	0.3864
	XGB	0.6818	0.6818	0.8636	0.5000	0.6333	0.3636
	MLP^*	0.6705	0.6705	0.6364	0.7045	0.6829	0.3409
Light	SVM	0.7841	0.7841	0.5909	0.9773	0.9630	0.5682
	NB	0.7273	0.7273	0.4773	0.9773	0.9545	0.4545
	KNN	0.8182	0.8182	0.7273	0.9091	0.8889	0.6364
	DT	0.7727	0.7727	0.6591	0.8864	0.8529	0.5455
	RF^†	0.7614	0.7614	0.8636	0.6591	0.7170	0.5227
	XGB	0.7273	0.7273	0.6364	0.8182	0.7778	0.4545
	MLP	0.7841	0.7841	0.5909	0.9773	0.9630	0.5682

Open in a new tab

The values for each performance metric were based on the feature set that achieved the highest AUROC in HV, as shown in Table 4. The numbers highlighted in bold typeface indicate the highest value in each column per light condition. The matricization strategy per classifier here is the highest shown in Table 4.

For MLP in light-on dataset, the full feature set of the $L 2 Norm$ index was used because it achieved the highest AUROC in HV.

^†

For RF in light-off dataset, the full feature set of $L 2 Norm$ index was used because it achieved the highest AUROC in HV.

In the light-on dataset, the classifier with the highest AUROC was NB (0.7386; Table 5). It also achieved the highest values in accuracy (0.7386), specificity (0.7500), precision (0.7442), and Kappa (0.4773). In addition, RF achieved the highest value in sensitivity (0.9545). In the light-off dataset, the classifier with the highest AUROC was KNN (0.8182). It achieved the highest values in accuracy (0.8182) and Kappa (0.6364). In addition, RF also achieved the highest sensitivity (0.8636). SVM, NB, and MLP achieved the highest specificity (0.9773). SVM and MLP also achieved the highest precision (0.9630). Therefore, the performance evaluated by other metrics was mostly consistent with that evaluated by the AUROC in light-on dataset, while the performance evaluated by other metrics varied from that evaluated by the AUROC in the light-off dataset.

4.8. Performance comparison in classification between features derived from matricization and CP decomposition

Although features from CP decomposition are less interpretable compared to features from matricization, features from CP decomposition might still be valuable. For instance, recent studies have proposed strategies to utilize features from CP decomposition for k-means clustering [72] and classifications [73]. Specifically, in the classification application, classifiers were trained using the training set projected into decomposed subspaces, using rank scores as features. Testing data were projected to the same subspace to evaluate the classification performances [73]. Even though using rank scores as features may be less interpretable because the contribution from features in the original 3DT space remains elusive, they may achieve better performance than that achieved by features derived from matricized data. Hence, we followed the established strategy for CP decomposition and compared their classification performance with features derived from matricization. Specifically, we conducted CP decomposition on the training set after confounding factor removal. The optimal rank number was determined by identifying the “elbow” in the scree plot of reconstruction error versus rank number [74]. Following the CP decomposition, we conducted cross-validation to identify the optimal hyperparameters for each classifier, using the same data partition as in the cross-validation with matricization (Fig 3). Afterwards, we removed confounding factors from the testing data, then projected the data into the same subspace using the Moore-Penrose Pseudoinverse strategy, the same projection method used in [73].

We compared the classification performances among the full feature sets from Index Construction, Feature Concatenation, and CP decomposition. The results for CV and HV are shown in S5 and S6 Tables, respectively. In CV, most classifiers achieved the highest AUROCs using Concat9Meas of Feature Concatenation (S5 Table). In the light-on dataset, 5 of the 7 classifiers achieved the highest AUROC using full features from matricization. Among these 5 classifiers, RF (0.8063), XGB (0.6960), and MLP (0.6110) used the features from Concat9Meas; NB (0.6038) used the features from Concat3Dist; DT (0.6574) used the features from Total Distance index. The two exceptions were SVM (0.7066) and KNN (0.6610), which achieved the highest AUROCs using features from CP decomposition. In the light-off dataset, 5 of the 7 classifiers achieved the highest AUROC using features from matricization. Among these 5 classifiers, NB (0.6558) and MLP (0.6901) used the features from Concat3Dist; RF (0.7633) used the features from Concat9Meas; DT (0.7360) used the features from the Total Distance index; XGB (0.7100) used the features from the $L 2 Norm$ index. The exceptions were SVM (0.6981) and KNN (0.6566), which achieved the highest AUROC using features from CP decomposition.

In HV, for the light-on dataset, 6 of the 7 classifiers achieved the highest AUROC using features from matricization (S6 Table). Among these 6 classifiers, DT (0.5795) and RF (0.7045) achieved the highest AUROC using features from Concat3Dist; XGB (0.6591) achieved the highest AUROC using features from Concat9Meas; KNN (0.6591) and MLP (0.6705) achieved the highest AUROC using features from the $L 2 Norm$ index; NB (0.7273) achieved the highest AUROC using features from $L 1 Norm$ index. The only exception was SVM (0.6023), which achieved the highest AUROC using features from CP decomposition. In the light-off dataset, 4 of the 7 classifiers achieved the highest AUROC using features from matricization. Among these 4 classifiers, DT (0.7045) achieved the highest AUROC using features from Concat9Meas; MLP (0.7159) achieved the highest AUROC using features from the Total Distance index; NB (0.6591) achieved the highest AUROC using features from the $L 1 Norm$ index; RF (0.7614) achieved the highest AUROC using features from the $L 2 Norm$ index. There were three exceptions: SVM (0.7273), KNN (0.6364), and XGB (0.6932) achieved the highest AUROC using features from CP decomposition. Altogether, the CV and HV results suggest that most classifiers achieved higher AUROCs using features from matricization than from CP decomposition.

5. Discussions

Contemporary neurobehavior research often analyzes multiple behavioral measurements collected from multiple samples over time. This results in MDT data that pose challenges to classical multivariate analysis, as such methods are designed to handle 2DT but not MDT data. Even though MDT data can be analyzed by CP decomposition, the resulting decomposed factors can be hard to interpret. Moreover, these results cannot be directly used for multivariate analysis, which is the preferred analysis approach of many neurobiologists. Consequently, they often transform the collected MDT data to 2DT data by various matricization strategies to enable the use of multivariate analysis and to ensure the interpretability of the analysis results. Various matricization methods can be applied to transform MDT features into a simplified 2DT space. However, some 2DT features are non-informative, and others may lose useful information during the matricization of the original data. These 2DT features must be carefully analyzed and selected for downstream multivariate analysis to optimize the analysis performance and facilitate the interpretation of underlying behavior dynamics. In this study, we assessed different matricization and feature-selection approaches and evaluated their impacts on the performance of a multivariate analysis.

Our study utilized a 3DT behavior dataset that contained zebrafish VMR collected from the wildtype and visually-impaired Q344X larvae. The VMR measures several aspects of swimming activity under a sudden light change, including $d u r$ , $c t$ , and $d i s t$ , at different activity levels. We first demonstrated the utility of matricization in interpretation by using Index Construction to create a Total Distance index, a commonly used index to plot and visualize VMR data [23,75,76]. The resulting plots of our VMR dataset (Fig 4) show a distinctive difference between the VMR of wildtype and Q344X, indicative of their differences in visual performance. In contrast, CP decomposition only reveals some patterns in each rank of respective factors (S1 Fig), but they do not easily facilitate the interpretation of the difference between the wildtype and Q344X VMR. We then evaluated the impact on multivariate analysis performance by various matricization and feature-selection approaches. Specifically, we conducted a classification task to distinguish between wildtype and Q334X larvae using matricized data. Several classifiers were built, and CV and HV were used to evaluate the performance of different matricization strategies and the selected 2DT features. Our results show that most classifiers achieved the highest AUROC in CV, using features from Feature Concatenation and selected by the embedded method or the union operation (Table 3). Specifically, in the light-on dataset, 4 of the 7 classifiers used features selected from Concat9Meas, and the other 3 used features selected from Concat3Dist. In the light-off dataset, almost all classifiers used features selected from Concat3Dist to achieve the highest performance. The superior performance of Feature Concatenation was also consistently seen in HV (Table 4). Their consistently high performance suggests that this specific combination of matricization strategy and feature selection methods (i.e., feature concatenation with embedded method/union operation) generally provides the best discriminative power for the VMR data. These results indicate several points: First, concatenating multiple measurements performed better than constructing an index. Different behavioral measurements in the MDT dataset may provide unique information and should be included in the downstream analysis whenever appropriate. Second, while most classifiers achieved better performance using features selected from Concat3Dist in the light-off dataset, some classifiers achieved better performance using features from Concat9Meas in the light-on dataset, suggesting that behavioral measurements impact classification performance differently, depending on light condition and the classifier. Hence, the matricization strategy should be chosen and evaluated based on the light condition in the VMR data and the choice of classifier. Nonetheless, conventional classifiers, such as NB and KNN, achieved the highest values in most performance metrics (Table 5) and outperformed more contemporary classifiers. This finding contradicts results from benchmark studies on other types of tabular data [69,77], where tree-based classifiers, particularly RF, outperformed conventional classifiers. It is possible that conventional classifiers performed better when only a few important features were selected and used for classification in our VMR data. For example, during CV on the light-on dataset (Table 3), RF performed worse when using features selected from Concat9Meas (0.7067, highest among all feature-selection methods) than when using the full feature set (0.8063). Similar results were observed during HV on the light-on dataset (Table 4), where RF achieved better performance using the full feature set (0.7045) than using selected features (0.6477, highest among all feature-selection methods). Indeed, RF employs a bootstrapping mechanism and an ensemble of decision trees, which can reduce the impact of noise from irrelevant features when the feature number is large [78-81]. Our results suggest that conventional classifiers benefited from feature selection and may outperform tree-based classifiers. Therefore, they should be considered a primary approach for analyzing VMR data with similar structures. Moreover, we also test if the HV performance measured by AUROC can be generalized to other performance metrics since they may reflect performance differently [77]. The results showed that the performance evaluated in AUROC was largely generalized to other metrics in the light-on dataset but not in the light-off dataset. These results are largely the same as a previous empirical evaluation in other tabular data, where certain classifiers performed better using certain metrics [77]. Nevertheless, AUROC can be used to evaluate classification performance in the light-on VMR dataset, while several other metrics should be used to evaluate the classification performance in the light-off VMR dataset, depending on the choice of classifier. In addition, we also compared the performances between features from matricization and features from CP decomposition. The features derived from matricization yielded superior classification performance across most classifiers evaluated in CV and HV. These results indicate that matricization outperforms CP decomposition in dimensionality reduction of 3DT data, including our VMR dataset. Combined with its superior interpretability, matricization is likely more effective than CP decomposition for the analysis of 3DT behavioral data.

The selected 2DT features reveal unique behavior dynamics between wildtype and Q344X larvae during light onset and offset. In the light-on VMR data, the filter method selected $d i s t a n c e . 1$ from the Total Distance index (Fig 5A) and $L 2 N o r m . 1$ from the $L 2 Norm$ index (Fig 5C). It also selected $i d u r . 1$ from the Concat9Meas (Fig 6B). The embedded method consistently selected small and inactive features derived from the first second, such as $d i s t a n c e . 1$ from the Total Distance index (Fig 7A), which incorporates the contribution of $i n a d i s t . 1$ from Concat3Dist and Concat9Meas, as well as $s m l d i s t . 1$ from Concat3Dist (Fig 8A and 8B). It also selected several small and inactive features derived from the first second unrelated to $d i s t$ , including $s m l c t . 1$ , $s m l d u r . 1$ , $i n a d u r . 1$ , and $i n a c t . 1$ from Concat9Mes (Fig 8B). Only a limited number of features from the later seconds were selected by the embedded methods, such as $L 2 N o r m . 41$ from $L 1 Norm$ index (Fig 7B) and $i n a d i s t . 131$ from Concat3Dist (Fig 8A). These results indicate that the most distinctive differences in the light-on VMR of the Q344X mutant were characterized by an absence of behavioral activity and/or a substantial reduction in the low-speed acute response during the first second, compared to the wildtype. Future studies of zebrafish visual mutants by scotopic light-on VMR should focus on analyzing the first-second behavior to maximize the discovery of neural mechanisms underlying the behavioral differences. In the light-off VMR, the filter method selected many features within the first 7 seconds in both datasets, matricized by Index Construction (Figs 5D-5E) and Feature Concatenation (Figs 6C and 6D). The embedded methods selected largely the same features as the filter method and revealed many of these features were associated with small and large activities (Figs 7D-7F & 8C-8D). For example, the top 3 features in Concat3Dist are $l a r d i s t . 2$ , $s m l d i s t . 3$ , and $l a r d i s t . 3$ (Fig 8C), whereas the top 3 features in Concat9Meas are $l a r d i s t . 2$ , $l a r c t . 2$ , and $l a r d u r . 2$ (Fig 8D). These results suggest that, compared to the Q344X mutant, more wildtype larvae displayed an acute startle response characterized by greater displacement sustained over longer durations. This is corroborated by our classification results, in which most classifiers achieved their highest AUROC with these features in both CV (Table 3) and HV (Table 4). Together, these results suggest that 2DT features in the early seconds after the light offset provide sufficient information to differentiate the wildtype from Q344X VMR. This is reasonable given that the normal light-off VMR typically exhibits the highest amplitude in the initial seconds (Fig 4B). Thus, the selected features from these initial seconds must provide a superior discriminative power. Even though the light-off VMR also tends to have a sustained response that can provide additional information, it may not be necessary for comparable discriminatory analyses.

The distinctive differences between our scotopic light-on and light-off VMR indicate that these responses are likely driven by different neural circuitries. Even though these scotopic circuitries are not well characterized, more is known about the ones driving the photopic startle response. Under a light flash (light onset), zebrafish larvae display a C-start response [82,83], consisting of a bend of about 70–80 degrees at around 183 ms, driven by the Mauthner cells [82]. Under a dark flash (light offset), zebrafish larvae display an O-bend[82,83] of approximately 150 degrees at around 408 ms [82]. This O-bend is not driven by the Mauthner cells but instead requires inputs from ventromedially-located hindbrain spinal projection neurons, including RoV3, MiV1, and MiV2 [84]. It has also been recently revealed that a distributed circuitry might control the habituation of O-bend [85,86]. The VMR is also shown by others and us to have both ocular and extra-ocular inputs [54,87,88]. The former involves contributions from both cones [89] and rods [54]. Different cone types seem to control different components of the photopic VMR [89], probably through different neural circuitries, as outlined above. However, the only rod system in the retina must control the observed scotopic VMR and innervate the corresponding neural circuitries driven by the different cone types under photopic stimulation.

Despite the new analytical and biological insights gained from our analyses, our study has several limitations. First, a few selected metrics were used for the filter methods, and one mechanism, RF, was used for the embedded method to illustrate the utility of our approach. This selection might introduce potential bias in the feature selection. More metrics and mechanisms could be used to select features to optimize the analysis performance [47]. Second, our analyses focused on matricization along the behavior-measurement dimension. Matricization, however, can be done on any dimension. For example, the time dimension can be matricized by constructing time-sensitive features [90,91], detecting time-series motifs [32,92], or calculating the area under the curve (AUC) index for each behavioral measurement in the time series [93]. Such matricization may also enhance downstream multivariate analysis performance and should be evaluated in future studies. Third, we acknowledge that more advanced feature selection methods exist beyond those we have tested. For instance, a recently developed filter method uses a Conditional-Weight Joint Relevance (CWJR) method to better capture the conditional relationships among features and achieve higher performance than other methods. Future work could compare and contrast the performance of CWJR and other feature selection methods in our dataset. Fourth, our analysis approach may not be suitable for all behavioral studies, especially when the experimental perturbation is known to influence a single behavioral measurement among multiple possible measurements. In this case, a simple univariate analysis is likely adequate for reaching a biological conclusion. For instance, in a study of hyperactive behavior mediated by mutations in GABA_A α subunits, the authors evaluated the effect of mutations on a touch-evoked response in the larvae by quantifying the resulting number of C-bends, a well-characterized initial phase of the C-start response [94]. Even though the larval behavior was captured by a high-speed camera and MDT data could presumably have been extracted and analyzed, such data were not essential for identifying the mutations in specific GABA subunits responsible for observed hyperactivity. Nonetheless, even in this case, it may be beneficial to analyze additional behavioral measurements other than the C-bend using our approach to maximize the likelihood of detecting behavioral effects caused by the mutations.

6. Conclusions

In conclusion, our study demonstrated the utility of different matricization approaches and feature-selection methods in analyzing MDT behavior data. Our evaluation approach can likely apply to other types of data as long as the data structure follows the MDT format. For instance, many zebrafish studies used calcium two-photon imaging to dissect neural circuitries driving several stereotypical behaviors in response to various stimuli [83-86,95,96]. The resulting data will contain multiple dimensions, including multiple behavioral measurements, time points, larvae, neurons, and their fluorescence readouts. Hence, these studies may benefit from using matricization strategies and feature-selection methods to identify relationships between the neuron firing pattern and the displayed behavior dynamics. In addition, neural-behavior studies conducted with animals other than zebrafish may also benefit from our studies, as their data also tend to consist of dimensions representing different trials, neurons, and time points of neural recordings within a single animal [17,28]. In time-series imaging studies, such as longitudinal magnetic resonance imaging [97], the data may consist of dimensions representing two dimensions of an image, samples, and time points [98,99]. To analyze these MDT data, we recommend the following workflow: Step 1) Preprocess and normalize the data using the standard procedures in the field; Step 2) Transform the MDT data to lower-order tensors by matricizing the dimensions of interest with the preferred approaches; Step 3) Select features from the matricized data with different feature-selection methods; Step 4) Evaluate the performance of specific combinations of matricization approaches and feature-selection methods; and Step 5) Apply the best-performing selected features for downstream multivariate analysis. By incorporating more information into the analysis, we anticipate that our workflow will uncover novel insights into complex neural behavior.

Supplementary Material

NIHMS2109698-supplement-1.docx^{(23.7KB, docx)}

NIHMS2109698-supplement-2.docx^{(244.1KB, docx)}

NIHMS2109698-supplement-3.docx^{(21KB, docx)}

NIHMS2109698-supplement-4.docx^{(22.1KB, docx)}

NIHMS2109698-supplement-5.docx^{(50.5KB, docx)}

NIHMS2109698-supplement-6.docx^{(40.1KB, docx)}

NIHMS2109698-supplement-7.docx^{(22.5KB, docx)}

NIHMS2109698-supplement-8.docx^{(22.2KB, docx)}

Highlights.

Many contemporary animal behavior data follow a multi-dimensional tensor format.
Multi-dimensional tensor can be analyzed by matricization or tensor decomposition.
Matricization yields more interpretable features than tensor decomposition.
In matricization, feature selection can extract important behavioral features.
Selected features increase classification performance across multiple classifiers.

Acknowledgments

We thank Mengrui Zhang for his thoughtful discussions.

Funding

This publication was made possible, in part, with support from the Indiana Clinical and Translational Sciences Institute to YFL and BW, funded in part by Grant Number UM1TR004402 from the National Institutes of Health, National Center for Advancing Translational Sciences, Clinical and Translational Sciences Award. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. JC, LF, and PM were partially supported by National Science Foundation grants DMS-2124493, DMS-2311297, DMS-2319279, DMS-2318809, and National Institutes of Health grants R01GM152814. MT was supported by Japan Society for the Promotion of Science grants JSPS KAKENHI 23K21480 and JSPS KAKENHI 24K22167.

Footnotes

Ethic Statement

All animal protocols were approved by the Purdue University Institutional Animal Care and Use Committee (Ref. 1201000592, Date: 03/07/2024).

The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Data Statement

The code supporting the results of this article is available on GitHub: <https://github.com/wang4537/VMR-Feature-Selection>. The dataset supporting the results of this article is available in the Harvard Dataverse repository: <https://doi.org/10.7910/DVN/B8HBU9>.

References

[1].Uddin LQ, Cognitive and behavioural flexibility: neural mechanisms and clinical considerations, Nat. Rev. Neurosci 22 (2021) 167–179. 10.1038/s41583-021-00428-w. [DOI] [PMC free article] [PubMed] [Google Scholar]
[2].Chen P, Hong W, Neural Circuit Mechanisms of Social Behavior, Neuron 98 (2018) 16–30. 10.1016/j.neuron.2018.02.026. [DOI] [PMC free article] [PubMed] [Google Scholar]
[3].Mench J, Why It Is Important to Understand Animal Behavior, ILAR J. 39 (1998) 20–26. 10.1093/ilar.39.1.20. [DOI] [PubMed] [Google Scholar]
[4].Gadagkar R, Observational study of animal behaviour: From instinct to intelligence, Curr. Sci 68 (1995) 185–196. [Google Scholar]
[5].Cuthill I, Field experiments in animal behaviour: methods and ethics, Anim. Behav 42 (1991) 1007–1014. 10.1016/S0003-3472(05)80153-8. [DOI] [Google Scholar]
[6].Warburton DM, Commentary on: “comprehensive observational assessment: Ia. A systematic, quantitative procedure for assessing the behavioral and physiologic state of the mouse.” Psychopharmacologia (1968) 13: 222-257, Psychopharmacology (Berl.) 163 (2002) 4–8. 10.1007/s00213-002-1119-x. [DOI] [PubMed] [Google Scholar]
[7].Bialek W, On the dimensionality of behavior, Proc. Natl. Acad. Sci 119 (2022) e2021860119. 10.1073/pnas.2021860119. [DOI] [PMC free article] [PubMed] [Google Scholar]
[8].Liu Y, Carmer R, Zhang G, Venkatraman P, Brown SA, Pang C-P, Zhang M, Ma P, Leung YF, Statistical Analysis of Zebrafish Locomotor Response, PLOS ONE 10 (2015) e0139521. 10.1371/journal.pone.0139521. [DOI] [PMC free article] [PubMed] [Google Scholar]
[9].Anderson MJ, Permutational Multivariate Analysis of Variance (PERMANOVA), in: Wiley StatsRef Stat. Ref. Online, John Wiley & Sons, Ltd, 2017: pp. 1–15. 10.1002/9781118445112.stat07841. [DOI] [Google Scholar]
[10].Liu Y, Ma P, Cassidy PA, Carmer R, Zhang G, Venkatraman P, Brown SA, Pang CP, Zhong W, Zhang M, Leung YF, Statistical Analysis of Zebrafish Locomotor Behaviour by Generalized Linear Mixed Models, Sci. Rep 7 (2017) 2937. 10.1038/s41598-017-02822-w. [DOI] [PMC free article] [PubMed] [Google Scholar]
[11].Kokel D, Bryan J, Laggner C, White R, Cheung CYJ, Mateus R, Healey D, Kim S, Werdich AA, Haggarty SJ, MacRae CA, Shoichet B, Peterson RT, Rapid behavior-based identification of neuroactive small molecules in the zebrafish, Nat. Chem. Biol 6 (2010) 231–237. 10.1038/nchembio.307. [DOI] [PMC free article] [PubMed] [Google Scholar]
[12].Hoffman EJ, Turner KJ, Fernandez JM, Cifuentes D, Ghosh M, Ijaz S, Jain RA, Kubo F, Bill BR, Baier H, Granato M, Barresi MJF, Wilson SW, Rihel J, State MW, Giraldez AJ, Estrogens Suppress a Behavioral Phenotype in Zebrafish Mutants of the Autism Risk Gene, CNTNAP2, Neuron 89 (2016) 725–733. 10.1016/j.neuron.2015.12.039. [DOI] [PMC free article] [PubMed] [Google Scholar]
[13].Kolda TG, Bader BW, Tensor Decompositions and Applications, SIAM Rev. 51 (2009) 455–500. 10.1137/07070111X. [DOI] [Google Scholar]
[14].De Lathauwer L, De Moor B, Vandewalle J, A Multilinear Singular Value Decomposition, SIAM J. Matrix Anal. Appl 21 (2000) 1253–1278. 10.1137/S0895479896305696. [DOI] [Google Scholar]
[15].Chen R, Yang D, Zhang C-H, Factor Models for High-Dimensional Tensor Time Series, J. Am. Stat. Assoc 117 (2022) 94–116. 10.1080/01621459.2021.1912757. [DOI] [Google Scholar]
[16].Carroll JD, Chang J-J, Analysis of individual differences in multidimensional scaling via an n-way generalization of “Eckart-Young” decomposition, Psychometrika 35 (1970) 283–319. 10.1007/BF02310791. [DOI] [Google Scholar]
[17].Williams AH, Kim TH, Wang F, Vyas S, Ryu SI, Shenoy KV, Schnitzer M, Kolda TG, Ganguli S, Unsupervised Discovery of Demixed, Low-Dimensional Neural Dynamics across Multiple Timescales through Tensor Component Analysis, Neuron 98 (2018) 1099–1115.e8. 10.1016/j.neuron.2018.05.015. [DOI] [PMC free article] [PubMed] [Google Scholar]
[18].Bro R, Smilde AK, Principal component analysis, Anal. METHODS 6 (2014) 2812–2831. 10.1039/c3ay41907j. [DOI] [Google Scholar]
[19].Yao F, Coquery J, Lê Cao K-A, Independent Principal Component Analysis for biologically meaningful dimension reduction of large biological data sets, BMC Bioinformatics 13 (2012) 24. 10.1186/1471-2105-13-24. [DOI] [PMC free article] [PubMed] [Google Scholar]
[20].Briscik M, Dillies M-A, Déjean S, Improvement of variables interpretability in kernel PCA, BMC Bioinformatics 24 (2023) 282. 10.1186/s12859-023-05404-y. [DOI] [PMC free article] [PubMed] [Google Scholar]
[21].Rabanser S, Shchur O, Günnemann S, Introduction to Tensor Decompositions and their Applications in Machine Learning, (2017). 10.48550/arXiv.1711.10781. [DOI] [Google Scholar]
[22].Zhang H, Zhang J, Wang R, Zhang Q, Spincemaille P, Nguyen TD, Wang Y, Efficient Folded Attention for 3D Medical Image Reconstruction and Segmentation, Proc. AAAI Conf. Artif. Intell 35 (2021) 10868–10876. 10.1609/aaai.v35i12.17298. [DOI] [Google Scholar]
[23].Ganzen L, Ko MJ, Zhang M, Xie R, Chen Y, Zhang L, James R, Mumm J, van Rijn RM, Zhong W, Pang CP, Zhang M, Tsujikawa M, Leung YF, Drug screening with zebrafish visual behavior identifies carvedilol as a potential treatment for an autosomal dominant form of retinitis pigmentosa, Sci. Rep 11 (2021) 11432. 10.1038/s41598-021-89482-z. [DOI] [PMC free article] [PubMed] [Google Scholar]
[24].Kroll F, Powell GT, Ghosh M, Gestri G, Antinucci P, Hearn TJ, Tunbak H, Lim S, Dennis HW, Fernandez JM, Whitmore D, Dreosti E, Wilson SW, Hoffman EJ, Rihel J, A simple and effective F0 knockout method for rapid screening of behaviour and other complex phenotypes, eLife 10 (2021) e59683. 10.7554/eLife.59683. [DOI] [PMC free article] [PubMed] [Google Scholar]
[25].Mosser EA, Chiu CN, Tamai TK, Hirota T, Li S, Hui M, Wang A, Singh C, Giovanni A, Kay SA, Prober DA, Identification of pathways that regulate circadian rhythms using a larval zebrafish small molecule screen, Sci. Rep 9 (2019) 12405. 10.1038/s41598-019-48914-7. [DOI] [PMC free article] [PubMed] [Google Scholar]
[26].Gustafsson FK, Danelljan M, Schön TB, Accurate 3D Object Detection using Energy-Based Models, (2023). 10.48550/arXiv.2012.04634. [DOI] [Google Scholar]
[27].Kiers HAL, Towards a standardized notation and terminology in multiway analysis, J. Chemom 14 (2000) 105–122. 10.1002/1099-128X(200005/06)14:3<105::AID-CEM582>3.0.CO;2-I. [DOI] [Google Scholar]
[28].Yu BM, Cunningham JP, Santhanam G, Ryu SI, Shenoy KV, Sahani M, Gaussian-Process Factor Analysis for Low-Dimensional Single-Trial Analysis of Neural Population Activity, J. Neurophysiol 102 (2009) 614–635. 10.1152/jn.90941.2008. [DOI] [PMC free article] [PubMed] [Google Scholar]
[29].Walker DM, Carmeli C, Pérez-Barbería FJ, Small M, Pérez-Fernández E, Inferring networks from multivariate symbolic time series to unravel behavioural interactions among animals, Anim. Behav 79 (2010) 351–359. 10.1016/j.anbehav.2009.11.011. [DOI] [Google Scholar]
[30].Yu Y, Becker T, Trinh LM, Behrisch M, SAXRegEx: Multivariate time series pattern search with symbolic representation, regular expression, and query expansion, Comput. Graph 112 (2023) 13–21. 10.1016/j.cag.2023.03.002. [DOI] [Google Scholar]
[31].Berman GJ, Choi DM, Bialek W, Shaevitz JW, Mapping the stereotyped behaviour of freely moving fruit flies, J. R. Soc. Interface 11 (2014) 20140672. 10.1098/rsif.2014.0672. [DOI] [PMC free article] [PubMed] [Google Scholar]
[32].Ghosh M, Rihel J, Hierarchical Compression Reveals Sub-Second to Day-Long Structure in Larval Zebrafish Behavior, eNeuro 7 (2020). 10.1523/ENEURO.0408-19.2020. [DOI] [PMC free article] [PubMed] [Google Scholar]
[33].Pudjihartono N, Fadason T, Kempa-Liehr AW, O’Sullivan JM, A Review of Feature Selection Methods for Machine Learning-Based Disease Risk Prediction, Front. Bioinforma 2 (2022) 927312. 10.3389/fbinf.2022.927312. [DOI] [PMC free article] [PubMed] [Google Scholar]
[34].Saeys Y, Inza I, Larrañaga P, A review of feature selection techniques in bioinformatics, Bioinformatics 23 (2007) 2507–2517. 10.1093/bioinformatics/btm344. [DOI] [PubMed] [Google Scholar]
[35].Dash M, Liu H, Feature selection for classification, Intell. Data Anal 1 (1997) 131–156. 10.1016/S1088-467X(97)00008-5. [DOI] [Google Scholar]
[36].Lal TN, Chapelle O, Weston J, Elisseeff A, Embedded Methods, in: Guyon I, Nikravesh M, Gunn S, Zadeh LA (Eds.), Feature Extr. Found. Appl, Springer, Berlin, Heidelberg, 2006: pp. 137–165. 10.1007/978-3-540-35488-8_6. [DOI] [Google Scholar]
[37].Jafari P, Azuaje F, An assessment of recently published gene expression data analyses: reporting experimental design and statistical factors, BMC Med. Inform. Decis. Mak 6 (2006) 27. 10.1186/1472-6947-6-27. [DOI] [PMC free article] [PubMed] [Google Scholar]
[38].Das S, Rai SN, Statistical Approach for Biologically Relevant Gene Selection from High-Throughput Gene Expression Data, Entropy 22 (2020) 1205. 10.3390/e22111205. [DOI] [PMC free article] [PubMed] [Google Scholar]
[39].Plackett RL, Karl Pearson and the Chi-Squared Test, Int. Stat. Rev. Rev. Int. Stat 51 (1983) 59–72. 10.2307/1402731. [DOI] [Google Scholar]
[40].Li J, Cheng K, Wang S, Morstatter F, Trevino RP, Tang J, Liu H, Feature Selection: A Data Perspective, ACM Comput. Surv 50 (2017) 94:1–94:45. 10.1145/3136625. [DOI] [Google Scholar]
[41].Alhaj TA, Siraj MM, Zainal A, Elshoush HT, Elhaj F, Feature Selection Using Information Gain for Improved Structural-Based Alert Correlation, PloS One 11 (2016) e0166017. 10.1371/journal.pone.0166017. [DOI] [PMC free article] [PubMed] [Google Scholar]
[42].Tibshirani R, Regression Shrinkage and Selection via the Lasso, J. R. Stat. Soc. Ser. B Methodol 58 (1996) 267–288. [Google Scholar]
[43].Breiman L, Random Forests, Mach. Learn 45 (2001) 5–32. 10.1023/A:1010933404324. [DOI] [Google Scholar]
[44].Quinlan JR, Induction of decision trees, Mach. Learn 1 (1986) 81–106. 10.1007/BF00116251. [DOI] [Google Scholar]
[45].Kudo M, Sklansky J, Comparison of algorithms that select features for pattern classifiers, Pattern Recognit. 33 (2000) 25–41. 10.1016/S0031-3203(99)00041-2. [DOI] [Google Scholar]
[46].Abu Khurma R, Aljarah I, Sharieh A, Abd Elaziz M, Damaševičius R, Krilavičius T, A Review of the Modification Strategies of the Nature Inspired Algorithms for Feature Selection Problem, Mathematics 10 (2022) 464. 10.3390/math10030464. [DOI] [Google Scholar]
[47].Bashir S, Khattak IU, Khan A, Khan FH, Gani A, Shiraz M, A Novel Feature Selection Method for Classification of Medical Data Using Filters, Wrappers, and Embedded Approaches, Complexity 2022 (2022) e8190814. 10.1155/2022/8190814. [DOI] [Google Scholar]
[48].Chen C-W, Tsai Y-H, Chang F-R, Lin W-C, Ensemble feature selection in medical datasets: Combining filter, wrapper, and embedded feature selection results, Expert Syst. 37 (2020) e12553. 10.1111/exsy.12553. [DOI] [Google Scholar]
[49].Nakao T, Tsujikawa M, Notomi S, Ikeda Y, Nishida K, The Role of Mislocalized Phototransduction in Photoreceptor Cell Death of Retinitis Pigmentosa, PLoS ONE 7 (2012) e32472. 10.1371/journal.pone.0032472. [DOI] [PMC free article] [PubMed] [Google Scholar]
[50].Westerfield M, The zebrafish book: a guide for the laboratory use of zebrafish (Danio rerio), 5th edition, Printed by the University of Oregon Press : Distributed by the Zebrafish International Resource Center, Eugene, OR, 2007. [Google Scholar]
[51].Hensley MR, Leung YF, A convenient dry feed for raising zebrafish larvae, Zebrafish 7 (2010) 219–231. 10.1089/zeb.2010.0652. [DOI] [PubMed] [Google Scholar]
[52].Emran F, Rihel J, Dowling JE, A behavioral assay to measure responsiveness of zebrafish to changes in light intensities, J. Vis. Exp. JoVE (2008) 923. 10.3791/923. [DOI] [PMC free article] [PubMed] [Google Scholar]
[53].Ganzen L, Venkatraman P, Pang CP, Leung YF, Zhang M, Utilizing Zebrafish Visual Behaviors in Drug Screening for Retinal Degeneration, Int. J. Mol. Sci 18 (2017) 1185. 10.3390/ijms18061185. [DOI] [PMC free article] [PubMed] [Google Scholar]
[54].Venkatraman P, Mills-Henry I, Padmanabhan KR, Pascuzzi P, Hassan M, Zhang J, Zhang X, Ma P, Pang CP, Dowling JE, Zhang M, Leung YF, Rods Contribute to Visual Behavior in Larval Zebrafish, Invest. Ophthalmol. Vis. Sci 61 (2020) 11. 10.1167/iovs.61.12.11. [DOI] [PMC free article] [PubMed] [Google Scholar]
[55].Xie R, Zhang M, Venkatraman P, Zhang X, Zhang G, Carmer R, Kantola SA, Pang CP, Ma P, Zhang M, Zhong W, Leung YF, Normalization of large-scale behavioural data collected from zebrafish, PloS One 14 (2019) e0212234. 10.1371/journal.pone.0212234. [DOI] [PMC free article] [PubMed] [Google Scholar]
[56].López Patiño MA, Yu L, Yamamoto BK, Zhdanova IV, Gender differences in zebrafish responses to cocaine withdrawal, Physiol. Behav 95 (2008) 36–47. 10.1016/j.physbeh.2008.03.021. [DOI] [PMC free article] [PubMed] [Google Scholar]
[57].Maleski ALA, Rosa JGS, Bernardo JTG, Astray RM, Walker CIB, Lopes-Ferreira M, Lima C, Recapitulation of Retinal Damage in Zebrafish Larvae Infected with Zika Virus, Cells 11 (2022) 1457. 10.3390/cells11091457. [DOI] [PMC free article] [PubMed] [Google Scholar]
[58].Gence L, Fernezelian D, Bringart M, Veeren B, Christophe A, Brion F, Meilhac O, Bascands J-L, Diotel N, Hypericum lanceolatum Lam. Medicinal Plant: Potential Toxicity and Therapeutic Effects Based on a Zebrafish Model, Front. Pharmacol 13 (2022). 10.3389/fphar.2022.832928. [DOI] [PMC free article] [PubMed] [Google Scholar]
[59].Li Y, Woodruff D, Yasuda T, Exponentially Improved Dimensionality Reduction for l1: Subspace Embeddings and Independence Testing, in: Proc. Thirty Fourth Conf. Learn. Theory, PMLR, 2021: pp. 3111–3195. https://proceedings.mlr.press/v134/li21c.html (accessed June 26, 2024). [Google Scholar]
[60].Guo H, Li Y, Mensah GK, Xu Y, Chen J, Xiang J, Chen D, Resting-State Functional Network Scale Effects and Statistical Significance-Based Feature Selection in Machine Learning Classification, Comput. Math. Methods Med 2019 (2019) 9108108. 10.1155/2019/9108108. [DOI] [PMC free article] [PubMed] [Google Scholar]
[61].Fay MP, Proschan MA, Wilcoxon-Mann-Whitney or t-test? On assumptions for hypothesis tests and multiple interpretations of decision rules, Stat. Surv 4 (2010) 1–39. 10.1214/09-SS051. [DOI] [PMC free article] [PubMed] [Google Scholar]
[62].Jordan R, Patel S, Hu H, Lyons-Weiler J, Efficiency Analysis of Competing Tests for Finding Differentially Expressed Genes in Lung Adenocarcinoma, Cancer Inform. 6 (2008) CIN.S791. 10.4137/CIN.S791. [DOI] [PMC free article] [PubMed] [Google Scholar]
[63].Holm S, A Simple Sequentially Rejective Multiple Test Procedure, Scand. J. Stat 6 (1979) 65–70. [Google Scholar]
[64].Breiman L, Friedman J, Olshen RA, Stone CJ, Classification and Regression Trees, Chapman and Hall/CRC, New York, 2017. 10.1201/9781315139470. [DOI] [Google Scholar]
[65].Sandri M, Zuccolotto P, Variable Selection Using Random Forests, in: Zani S, Cerioli A, Riani M, Vichi M (Eds.), Data Anal. Classif. Forw. Search, Springer, Berlin, Heidelberg, 2006: pp. 263–270. 10.1007/3-540-35978-8_30. [DOI] [Google Scholar]
[66].Nembrini S, König IR, Wright MN, The revival of the Gini importance?, Bioinformatics 34 (2018) 3711–3718. 10.1093/bioinformatics/bty373. [DOI] [PMC free article] [PubMed] [Google Scholar]
[67].Gao Y, Zhang G, Jelfs B, Carmer R, Venkatraman P, Ghadami M, Brown SA, Pang CP, Leung YF, Chan RHM, Zhang M, Computational classification of different wild-type zebrafish strains based on their variation in light-induced locomotor response, Comput. Biol. Med 69 (2016) 1–9. 10.1016/j.compbiomed.2015.11.012. [DOI] [PubMed] [Google Scholar]
[68].Geurts P, Irrthum A, Wehenkel L, Supervised learning with decision tree-based methods in computational and systems biology, Mol. Biosyst 5 (2009) 1593. 10.1039/b907946g. [DOI] [PubMed] [Google Scholar]
[69].Grinsztajn L, Oyallon E, Varoquaux G, Why do tree-based models still outperform deep learning on typical tabular data?, Adv. Neural Inf. Process. Syst 35 (2022) 507–520. [Google Scholar]
[70].Pang B, Nijkamp E, Wu YN, Deep Learning With TensorFlow: A Review, J. Educ. Behav. Stat 45 (2020) 227–248. 10.3102/1076998619872761. [DOI] [Google Scholar]
[71].Sheela KG, Deepa SN, Review on Methods to Fix Number of Hidden Neurons in Neural Networks, Math. Probl. Eng 2013 (2013) 425740. 10.1155/2013/425740. [DOI] [Google Scholar]
[72].Qi J, Zhang J, Zhang Y, Jiang X, Cai Z, Tensor Decomposition Based Latent Feature Clustering for Hyperspectral Band Selection, in: ICASSP 2023 - 2023 IEEE Int. Conf. Acoust. Speech Signal Process. ICASSP, 2023: pp. 1–5. 10.1109/ICASSP49357.2023.10096731. [DOI] [Google Scholar]
[73].Amin MR, Hasan M, Arnab SP, DeGiorgio M, Tensor Decomposition-based Feature Extraction and Classification to Detect Natural Selection from Genomic Data, Mol. Biol. Evol 40 (2023) msad216. 10.1093/molbev/msad216. [DOI] [PMC free article] [PubMed] [Google Scholar]
[74].Heng Q, Chi EC, Liu Y, Robust Low-rank Tensor Decomposition with the $L 2$ Criterion, Technometrics J. Stat. Phys. Chem. Eng. Sci 65 (2023) 537–552. 10.1080/00401706.2023.2200541. [DOI] [PMC free article] [PubMed] [Google Scholar]
[75].Beppi C, Beringer G, Straumann D, Bögli SY, Light-stimulus intensity modulates startle reflex habituation in larval zebrafish, Sci. Rep 11 (2021) 22410. 10.1038/s41598-021-00535-9. [DOI] [PMC free article] [PubMed] [Google Scholar]
[76].Beppi C, Straumann D, Bögli SY, A model-based quantification of startle reflex habituation in larval zebrafish, Sci. Rep 11 (2021) 846. 10.1038/s41598-020-79923-6. [DOI] [PMC free article] [PubMed] [Google Scholar]
[77].Caruana R, Niculescu-Mizil A, An empirical comparison of supervised learning algorithms, in: Proc. 23rd Int. Conf. Mach. Learn., Association for Computing Machinery, New York, NY, USA, 2006: pp. 161–168. 10.1145/1143844.1143865. [DOI] [Google Scholar]
[78].Li S, Harner EJ, Adjeroh D, Random KNN feature selection - a fast and stable alternative to Random Forests, BMC Bioinformatics 12 (2011). 10.1186/1471-2105-12-450. [DOI] [PMC free article] [PubMed] [Google Scholar]
[79].Yao D, Zhan X, Zhan X, Kwoh CK, Peng L, Wang J, A random forest based computational model for predicting novel lncRNA-disease associations, BMC Bioinformatics 21 (2020). 10.1186/s12859-020-3458-1. [DOI] [PMC free article] [PubMed] [Google Scholar]
[80].Tandra S, Manashty A, Probabilistic Feature Selection for Interpretable Random Forest Model, Adv. Intell. Syst. Comput (2021). 10.1007/978-3-030-73103-8_50. [DOI] [Google Scholar]
[81].Probst P, Wright M, Boulesteix A-L, Hyperparameters and tuning strategies for random forest, WIREs Data Min. Knowl. Discov 9 (2019). 10.1002/widm.1301. [DOI] [Google Scholar]
[82].Burgess HA, Granato M, Modulation of locomotor activity in larval zebrafish during light adaptation, J. Exp. Biol 210 (2007) 2526–2539. 10.1242/jeb.003939. [DOI] [PubMed] [Google Scholar]
[83].Marques JC, Lackner S, Félix R, Orger MB, Structure of the Zebrafish Locomotor Repertoire Revealed with Unsupervised Behavioral Clustering, Curr. Biol 28 (2018) 181–195.e5. 10.1016/j.cub.2017.12.002. [DOI] [PubMed] [Google Scholar]
[84].Huang K-H, Ahrens MB, Dunn TW, Engert F, Spinal Projection Neurons Control Turning Behaviors in Zebrafish, Curr. Biol 23 (2013) 1566–1573. 10.1016/j.cub.2013.06.044. [DOI] [PMC free article] [PubMed] [Google Scholar]
[85].Lamiré LA, Haesemeyer M, Engert F, Granato M, Randlett O, Functional and pharmacological analyses of visual habituation learning in larval zebrafish, eLife 12 (2023) RP84926. 10.7554/eLife.84926. [DOI] [PMC free article] [PubMed] [Google Scholar]
[86].Randlett O, Haesemeyer M, Forkin G, Shoenhard H, Schier AF, Engert F, Granato M, Distributed Plasticity Drives Visual Habituation Learning in Larval Zebrafish, Curr. Biol 29 (2019) 1337–1345.e4. 10.1016/j.cub.2019.02.039. [DOI] [PMC free article] [PubMed] [Google Scholar]
[87].Emran F, Rihel J, Adolph AR, Wong KY, Kraves S, Dowling JE, OFF ganglion cells cannot drive the optokinetic reflex in zebrafish, Proc. Natl. Acad. Sci. U. S. A 104 (2007) 19126–19131. 10.1073/pnas.0709337104. [DOI] [PMC free article] [PubMed] [Google Scholar]
[88].Fernandes AM, Fero K, Arrenberg AB, Bergeron SA, Driever W, Burgess HA, Deep Brain Photoreceptors Control Light-Seeking Behavior in Zebrafish Larvae, Curr. Biol 22 (2012) 2042–2047. 10.1016/j.cub.2012.08.016. [DOI] [PMC free article] [PubMed] [Google Scholar]
[89].Burton CE, Zhou Y, Bai Q, Burton EA, Spectral properties of the zebrafish visual motor response, Neurosci. Lett 646 (2017) 62–67. 10.1016/j.neulet.2017.03.002. [DOI] [PMC free article] [PubMed] [Google Scholar]
[90].Gao Y, Chan RHM, Chow TWS, Zhang L, Bonilla S, Pang C-P, Zhang M, Leung YF, A High-Throughput Zebrafish Screening Method for Visual Mutants by Light-Induced Locomotor Response, IEEE/ACM Trans. Comput. Biol. Bioinform 11 (2014) 693–701. 10.1109/TCBB.2014.2306829. [DOI] [PubMed] [Google Scholar]
[91].Christ M, Braun N, Neuffer J, Kempa-Liehr AW, Time Series FeatuRe Extraction on basis of Scalable Hypothesis tests (tsfresh – A Python package), Neurocomputing 307 (2018) 72–77. 10.1016/j.neucom.2018.03.067. [DOI] [Google Scholar]
[92].Mueen A, Time series motif discovery: dimensions and applications, WIREs Data Min. Knowl. Discov 4 (2014) 152–159. 10.1002/widm.1119. [DOI] [Google Scholar]
[93].Carbaugh CM, Widder MW, Phillips CS, Jackson DA, DiVito VT, Schalie WH, Glover KP, Assessment of zebrafish embryo photomotor response sensitivity and phase-specific patterns following acute- and long-duration exposure to neurotoxic chemicals and chemical weapon precursors, J. Appl. Toxicol 40 (2020) 1272–1283. 10.1002/jat.3984. [DOI] [PMC free article] [PubMed] [Google Scholar]
[94].Barnaby W, Dorman Barclay HE, Nagarkar A, Perkins M, Teicher G, Trapani JG, Downes GB, GABAA α subunit control of hyperactive behavior in developing zebrafish, Genetics 220 (2022) iyac011. 10.1093/genetics/iyac011. [DOI] [PMC free article] [PubMed] [Google Scholar]
[95].Ahrens MB, Huang K-H, Narayan S, Mensh BD, Engert F, Two-photon calcium imaging during fictive navigation in virtual environments, Front. Neural Circuits 7 (2013). 10.3389/fncir.2013.00104. [DOI] [PMC free article] [PubMed] [Google Scholar]
[96].Ahrens MB, Li JM, Orger MB, Robson DN, Schier AF, Engert F, Portugues R, Brain-wide neuronal dynamics during motor adaptation in zebrafish, Nature 485 (2012) 471–477. 10.1038/nature11057. [DOI] [PMC free article] [PubMed] [Google Scholar]
[97].Sen B, Parhi KK, Constrained Tensor Decomposition Optimization With Applications To Fmri Data Analysis, in: 2018 52nd Asilomar Conf. Signals Syst. Comput., 2018: pp. 1923–1928. 10.1109/ACSSC.2018.8645427. [DOI] [Google Scholar]
[98].Monti MM, Statistical Analysis of fMRI Time-Series: A Critical Review of the GLM Approach, Front. Hum. Neurosci 5 (2011) 28. 10.3389/fnhum.2011.00028. [DOI] [PMC free article] [PubMed] [Google Scholar]
[99].Brand L, Nichols K, Wang H, Huang H, Shen L, for the Alzheimer’s Disease Neuroimaging Initiative, Predicting Longitudinal Outcomes of Alzheimer’s Disease via a Tensor-Based Joint Classification and Regression Model, in: Biocomput. 2020, WORLD SCIENTIFIC, Kohala Coast, Hawaii, USA, 2019: pp. 7–18. 10.1142/9789811215636_0002. [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

NIHMS2109698-supplement-1.docx^{(23.7KB, docx)}

NIHMS2109698-supplement-2.docx^{(244.1KB, docx)}

NIHMS2109698-supplement-3.docx^{(21KB, docx)}

NIHMS2109698-supplement-4.docx^{(22.1KB, docx)}

NIHMS2109698-supplement-5.docx^{(50.5KB, docx)}

NIHMS2109698-supplement-6.docx^{(40.1KB, docx)}

NIHMS2109698-supplement-7.docx^{(22.5KB, docx)}

NIHMS2109698-supplement-8.docx^{(22.2KB, docx)}

PERMALINK

Tensor analysis of animal behavior by matricization and feature selection

Beichen Wang

Jiazhang Cai

Luyang Fang

Motokazu Tsujikawa

Ping Ma

Yuk Fai Leung

Abstract

1. Introduction

Fig 1. Schematic structures of behavior data in multivariate and multidimensional tensor form.

2. Background

2.1. Multi-dimensional tensor

2.2. Feature selection

Table 1. Three methods for feature selection.

3. Methods

3.1. Zebrafish husbandry and breeding

3.2. Visual-motor response assay

3.3. VMR data processing and analysis

3.3.1. General

3.3.2. VMR data structure

Fig 2. Schematics of the 3DT from the zebrafish behavior data and matricization strategies.

3.3.3. Matricizing the 3DT VMR data

3.3.4. Feature selection

3.3.5. Evaluation of feature selection and matricization performance

Fig 3. Evaluation of different matricization and feature-selection methods.

4. Results

4.1. The 3DT VMR dataset

Fig 4. Wildtype and Q344X VMR illustrated by the Total Distance index.

4.2. Feature Interpretation

4.3. Feature selection—filter method

Fig 5. Feature selection by the filter method after Index Construction in the VMR datasets.

Fig 6. Feature selection by the filter method after Feature Concatenation in the VMR datasets.

4.4. Feature selection—embedded method

Fig 7. Feature selection by the embedded method after Index Construction in the VMR datasets.

Fig 8. Feature selection by the embedded method after Feature Concatenation in the VMR datasets.

4.5. Aggregation of selected 2DT features and performance evaluation

Table 2. The number of features selected by different approaches.

Table 3.

4.6. Generalization of the matricization approach and selected feature set to new data

Table 4. The highest AUROC of HV in Index Construction and Feature Concatenation.

4.7. Performance evaluation of the matricization approach and selected features by other performance metrics

Table 5.

4.8. Performance comparison in classification between features derived from matricization and CP decomposition

5. Discussions

6. Conclusions

Supplementary Material

Highlights.

Acknowledgments

Funding

Footnotes

Data Statement

References

Associated Data

Supplementary Materials

Data Availability Statement

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases