anomalyDetection: Implementation of Augmented Network Log Anomaly Detection Procedures

As the number of cyber-attacks continues to grow on a daily basis, so does the delay in threat detection. For instance, in 2015, the Office of Personnel Management discovered that approximately 21.5 million individual records of Federal employees and contractors had been stolen. On average, the time between an attack and its discovery is more than 200 days. In the case of the OPM breach, the attack had been going on for almost a year. Currently, cyber analysts inspect numerous potential incidents on a daily basis, but have neither the time nor the resources available to perform such a task. anomalyDetection aims to curtail the time frame in which anomalous cyber activities go unnoticed and to aid in the efficient discovery of these anomalous transactions among the millions of daily logged events by i) providing an efficient means for pre-processing and aggregating cyber data for analysis by employing a tabular vector transformation and handling multicollinearity concerns; ii) offering numerous built-in multivariate statistical functions such as Mahalanobis distance, factor analysis, principal components analysis to identify anomalous activity, iii) incorporating the pipe operator (%>%) to allow it to work well in the tidyverse workflow. Combined, anomalyDetection offers cyber analysts an efficient and simplified approach to break up network events into time-segment blocks and identify periods associated with suspected anomalies for further evaluation.


Introduction
Organizations worldwide rely heavily on the systems of cyberspace, and the wider Internet as a whole, for commerce, defense operations, infrastructure control systems, financial management, transportation, and other critical services. Unfortunately, the number of cyber-attacks are growing on a daily basis and the ability of organizations to spot anomalous cyber activitiy is becoming more and more delayed . On average, the time between an attack and its discovery is more than 200 days (Koerner, 2016). In the case of the 2015 Office of Personnel Management breach, in which approximately 21.5 million individual records of Federal employees and contractors had been stolen, the attack had been going on for almost a year prior to the anomalous activity being identified (U. S. Office of Personnel Management, 2015).
Cyber analysts inspect numerous potential incidents on a daily basis, but lack the time and resources to scour the high volumes of data in an efficient manner to identify anomalous activity worth further investigation (Samuelson, 2016). Firewalls and intrusion detection and prevention systems (IDPS) are one line of defense in identifying and stopping suspicious internet traffic. When a suspicious event occurs, these devices generate a log file containing details of what preprogrammed rules were violated and how it was handled (Goodall et al., 2009). Such log files contain details of the event, (i.e. source and destination IP addresses, port numbers, and protocols), but not the packet and data that led to the event. A primary activity of cyber analysts is the analysis of these log files to detect anomalies . Although a reduced form of anomalous cyber data, these data sets can still represent millions of cyberspace transactions per minute (Jayathilake, 2012). Unfortunately, analysis of these log files has, historically, been heavily manual in nature and often leverages subject matter expertise to find possible threats in logged events to further investigate (Goodall et al., 2009;Jayathilake, 2012;Samuelson, 2016). This approach is inefficient and can suffer from biased, subjective assessments (Zamani, 2013). Moreover, much of the related research has focused on anomaly detection at the device/software level (i.e. Lazarevic et al., 2003;Denning, 1987;Garcia-Teodoro et al., 2009), with little exploration into anomaly detection in the log files generated from the preexisting devices or software (i.e McDonald et al., 2012;Winding et al., 2006;Breier and Branišová, 2015). Consequently, efficient analytic approaches are desirable to help detect anomalous activity in cyber network log data .
This research introduces the anomalyDetection package  to provide cyber analysts efficient means for performing anomaly detection in log files. The purpose of the package is to make identifying abnormal activity more efficient so that cyber analysts can spend more time researching the potential threat. It is important to note that there is no guarantee that anomalous activity is evidence of malicious cyber activity; however, identification of anomalous activity provides cyber security experts a starting point in their search for undetected malicious activity. anomalyDetection simplifies this process. This paper proceeds as follows. First, we introduce the anomalyDetection functions and the The R Journal Vol. 9/2, December 2017 ISSN 2073-4859 security_logs data set that will be leveraged for the illustrative examples. Next, we demonstrate how anomalyDetection can pre-process a data set for follow on analysis. This includes converting a data set with non-numeric attributes into numeric data using a tabulated state vector approach. anomalyDetection further prepares this tabulated state vector by inspecting and correcting for multicollinearity. We then illustrate how anomalyDetection provides efficient multivariate statistical analysis processes to help identify anomalous activity. Last, we end with some concluding remarks.
anomalyDetection functions anomalyDetection provides 13 functions to aid in the detection of potential cyber anomalies, which are listed in Table 1. The package also incorporates the pipe operator (%>%) from the magrittr package (Bache and Wickham, 2014) for streamlining function composition. To illustrate the functionality of anomalyDetection we will use the security_logs data that mimics common information that appears in security logs and comes with anomalyDetection. Note that we also load the tidyverse package (Wickham, 2017) for common manipulation and visualization tasks.

mahalanobis_distance()
Calculates the distance between the elements in data and the mean vector of the data for outlier detection. bd_row() Indicates which variables in data are driving the Mahalanobis distance for a specific row, relative to the mean vector of the data. horns_curve() Computes Horn's Parallel Analysis to determine the factors to retain within a factor analysis. factor_analysis() Reduces the structure of the data by relating the correlation between variables to a set of factors, using the eigen-decomposition of the correlation matrix. factor_analysis_results() Provides easy access to factor analysis results. kaisers_index() Computes scores designed to assess the quality of a factor analysis solution. It measures the tendency towards unifactoriality for both a given row and the entire matrix as a whole.

principal_components()
Relates the data to a set of a components through the eigendecomposition of the correlation matrix, where each component explains some variance of the data. principal_components_results() Provides easy access to principal component analysis results. get_all_factors() finds all factor pairs for a given integer.

Data pre-processing
In order to develop a statistical framework for firewall log analysis, data pre-processing is necessary prior to applying any multivariate analytic techniques. To assist in this process anomalyDetection offers two main approaches to pre-process log file data -aggregating the data into a tabulated state vector and managing multicollinearity concerns.
First, we can employ the tabulated vector approach introduced by . This approach transforms the security log data into unique counts of data attributes based on pre-defined time blocks. Therefore, as each time block is generated, the categorical fields are separated by their levels and a count of occurrences for each level are recorded into a vector. All numerical fields, such as bytes in and bytes out, are recorded as a summation within the time block. The result is what we call a state vector matrix.
Thus, for our security_logs data we can create our state vector matrix based on our data being divided into 10 time blocks. What results is the summary of instances for each categorical level in our data for each time block. Consequently, row one represents the first time block and there were 2 instances of CISCO as the device vendor, 1 instance of IBM, etc. By adjusting the block_length, level_limit, and level_keep arguments, the user can refine the level of aggregation and variables to retain and analyze. The state vector matrix provides us with a numerical construct to analyze our log file data; however, prior to proceeding with any multivariate statistical analyses we should inspect the state vector for multicollinearity, to avoid issues such as matrix singularity, rank deficiency, and strong correlation values, and remove any columns that pose an issue. We can use mc_adjust() to handle issues with multicollinearity by first removing any columns whose variance is close to or less than a minimum level of variance (min_var). Then, it removes linearly dependent columns. Finally, it removes any columns that have a high absolute correlation value equal to or greater than that defined by the user (max_cor).
(state_vec <-security_logs %>% tabulate_state_vector (10)  By default, mc_adjust() removes all columns that violate the variance, dependency, and correlation thresholds. Alternatively, we can use action = "select" as an argument, which provides interactivity where the user can select the variables that violate the correlation threshold that they would like to remove.

Multivariate statistical analyses
With our data adjusted for multicollinearity we can now proceed with multivariate analyses to identify anomalies in our log file. First we'll use the mahalanobis_distance() function to compare the distance between each observation by its distance from the data mean, independent of scale. This is computed as x is the mean vector of the data,x = (x 1 , . . . ,x p ), and C −1 is the inverse data covariance matrix. Here, we include output = "both" to return both the Mahalanobis distance and the absolute breakdown distances and normalize = TRUE so that we can compare relative magnitudes across our data. We can use this information in a modified heatmap visualization (Figure 1) to identify outlier values across our security log attributes and time blocks. Brighter columns represent time blocks deserving greater attention and further investigation. Larger circles represent variables within a time block that have more anomalous activity. Thus, the larger and brighter the dot the more significant the outlier is and the more it deserves attention.
state_vec %>% mahalanobis_distance("both", normalize = TRUE) %>% as_tibble %>% dplyr::mutate(Block = 1:n()) %>% gather(Variable, BD, -c(MD, Block)) %>% ggplot(aes(factor(Block), Variable, color = MD, size = BD)) + geom_point() We can build onto this with the bd_row() function to identify which security log attributes in the data are driving the Mahalanobis distance. bd_row() measures the relative contribution of each variable, x i , to MD by computing where C ii is the variance of x i . Furthermore, bd_row() will look at a specified row and rank-order the columns by those that are driving the Mahalanobis distance. For example, the plot above identified block 17 as having the largest Mahalanobis distance suggesting some abnormal activity may be occurring during that time block. We can drill down into that block and look at the top 10 security log attributes that are driving the Mahalanobis distance as these may be areas that require further investigation. Next, we can use factor analysis as a dimensionality reduction technique to identify the underlying structure of the data and identify factors (features) in the data that appear abnormal. Factor analysis relates the correlations between variables through a set of factors to link together seemingly unrelated variables. The basic factor analysis model is where X is the vector of responses X = (x 1 , . . . , x p ), f are the common factors f = ( f 1 , . . . , f q ), e is the unique factors e = (e 1 , . . . , e p ), and Λ is the factor loadings. Factor loadings are correlations between the factors and the original data and can thus range from -1 to 1, which indicate how much that factor affects each variable. Values close to 0 imply a weak effect on the variable. For the desired results, anomalyDetection uses the correlation matrix in its factor analysis computation.
A factor loadings matrix can be computed to understand how each original data variable is related to the resultant factors. This can be computed as where λ 1 is the eigenvalue for each factor, e i is the eigenvector for each factor, and p is the number of columns. Factor scores are used to examine the behavior of the observations relative to each factor and can be used to identify anomaly detection. Factor scores are calculated aŝ where X s is the standardized observations, R −1 is the inverse of the correlation matrix, and Λ is the factor loadings matrix. To simplify the results for interpretation, the factor loadings can undergo an orthogonal or oblique rotation. Orthogonal rotations assume independence between the factors while oblique rotations allow the factors to correlate. anomalyDetection utilizes the most common rotation option known as varimax. Varimax rotates the factors orthogonally to maximize the variance of the squared factor loadings which forces large factors to increase and small ones to decrease, providing easier interpretation.
To begin using factor analysis, the dimensions of the reduced state vector matrix are first passed to horns_curve(), which computes Horn's Parallel Analysis (Horn, 1965) to determine the factors to retain within a factor analysis. Next, the dimensionality is determined by finding the eigenvalues of the correlation matrix of the state vector matrix and retaining only those factors whose eigenvalues are greater than or equal to those produced by horns_curve(). We use factor_analysis() to reduce the state vector matrix into resultant factors. The factor_analysis() function generates a list containing five outputs: $fa_loadings numerical matrix with the original factor loadings $fa_scores numerical matrix with the row scores for each factor $fa_loadings_rotated numerical matrix with the varimax rotated factor loadings $fa_scores_rotated numerical matrix with the row scores for each varimax rotated factor $num_factors : numeric vector identifying the number of factors state_vec %>% horns_curve() %>% factor_analysis(state_vec, hc_points = .) %>% names() # [1] "fa_loadings" "fa_scores" "fa_loadings_rotated" # [4] "fa_scores_rotated" "num_factors" For easy access to these results we can use the factor_analysis_results() parsing function. The factor_analysis_results() function will parse the results either by their list name or by location. For instance to extract the rotated factor scores you can use factor_analysis_results(data, results = fa_scores_rotated) or factor_analysis_results(data, results = 4) as demonstrated below. To evaluate the quality of a factor analysis solution, Kaiser (Kaiser, 1974) proposed the Index of Factorial Simplicity (IFS). The IFS is computed as where q is the number of factors, j the row index, s the column index, and λ js is the value in the loadings matrix. Furthermore, Kaiser created the following evaluations of the score produced by the IFS as shown below:

In the .80s Meritorious
In the .70s Middling

In the .60s Mediocre
In the .50s Miserable Less than .50 : Unacceptable Thus, to assess the quality of our factor analysis results we apply kaisers_index() to the rotated factor loadings and, as the results show below, our output value of 0.702 suggests that our results are "middling". state_vec %>% horns_curve() %>% factor_analysis(data = state_vec, hc_points = .) %>% factor_analysis_results(fa_loadings_rotated) %>% kaisers_index() # [1] 0.7018006 Furthermore, Figure 2 visualizes the factor analysis results to show the correlation between the columns of the reduced state vector to the rotated factor loadings. Strong negative correlations are depicted as red while strong positive correlations are shown as blue. This helps to identify which factors are correlated with each security log data attribute. Furthermore, this helps to identify two or more security log data attributes that appear to have relationships with their occurrences. For example, this shows that Russia is highly correlated with IP address 223.70.128 since both these attributes are strongly correlated with factor 5. If there is an abnormally large amount of instances with Russian occurrences this would be the logical IP address to start investigating.
We can also visualize the rotated factor score plots as in Figure 3 to see which time blocks appear to be outliers and deserve closer attention.
An alternative, yet similar approach to factor analysis is principal component analysis. These two approaches can produce similar outcomes, especially when the error component in equation 4 is close to zero (Fabrigar et al., 1999). However, the results often differ and there are important distinctions in the interpretation of these results (Park et al., 2002). First, a primary difference between the two approaches is that factor analysis estimates errors while principal component analysis does not. This indicates that principal component analysis assumes that the measurement is without error. Second, the goal in factor analysis is to explain the covariances or correlations between the variables. Therefore anomaly detection using factor analysis will identify time blocks and variable attributes in which their abnormal behavior is highly correlated to one another. This allows you to identify latent features in the data. By contrast, the goal of principal component analysis is to explain as much of the total variance in the variables as possible. Therefore, if your goal is to reduce the log file variables into a composite component for further analysis, principal component analysis would be appropriate. To maintain clarity between these two approaches the following discussion leverages different notation.
The first principal component of a set of features X 1 , X 2 , . . . , X p is the normalized linear combination of the features that has the largest variance. By normalized , we mean that ∑ p j=1 φ 2 j1 = 1. We refer to the elements φ 11 , . . . , φ p1 as the loadings of the first principal component; together, the loadings make up the principal component loading vector, φ 1 = (φ 11 , φ 21 , . . . , φ p1 ) T . The loadings are constrained so that their sum of squares is equal to one, since otherwise setting these elements to be arbitrarily large in absolute value could result in an arbitrarily large variance. After the first principal component Z 1 of the features has been determined, we can find the second principal component Z 2 . The second principal component is the linear combination of X 1 , . . . , X p that has maximal variance out of all linear combinations that are uncorrelated with Z 1 . The second principal component scores z 12 , z 22 , . . . , z n2 take the form where φ 2 is the second principal loading vector, with elements φ 12 , φ 22 , . . . , φ p2 . This continues until all principal components have been computed. Therefore anomaly detection using PCA will maximize the difference in behaviors across time blocks and variable attributes. Thus, identifying anomalies with PCA will identify those attributes that behave very differently than all the other features. To perform a principal components analysis we use principal_components() which will create a list containing: $pca_sdev the standard deviations of the principal components (i.e., the square roots of the eigenvalues of the covariance/correlation matrix, though the calculation is actually done with the singular values of the data matrix). $pca_loadings the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). $pca_rotated the value of the rotated data (the centered, and scaled if requested, data multiplied by the rotation matrix) is returned. $pca_center the centering used. $pca_scale a logical response indicating whether scaling was used.
principal_components(state_vec) %>% names # [1] "pca_sdev" "pca_loadings" "pca_rotated" "pca_center" # [5] "pca_scale" For easy access to these results we can use the principal_components_result parsing function. The principal_components_result will parse the results either by their list name or by location. For example, to extract the computed component scores as outlined in Eq. 8 you can use We can then follow the principal components analysis with similar visualization activities as performed post-factor analysis to identify features that exhibit abnormal behavior. Since visualizing principal components analysis to identify anomolies mirrors that which we performed in the factor analysis section, we will leave this to the reader as an independent exercise.

Conclusion
Cyber attacks continue to be a growing concern for organizations. Unfortunately, the process of analyzing log files has, historically, been unorganized and lacked efficient approaches. The presented anomalyDetection package makes the log file analysis process more efficient and facilitates the identification and analysis of anomalies within log files.
First, anomalyDetection improves the pre-processing of cyber data. The package offers functions that help to narrow down abnormal behavior by aggregating internet traffic data into customizable time blocks. Aggregated activity at a higher-level time block should be easier to analyze while still offering a map to suspicious areas to drill down using smaller time blocks. For very large data sets, the analyst can tune the function parameters to start with less blocks and more aggregated data and then iteratively drill down into less aggregated data. Furthermore, anomalyDetection improves the process of adjusting for multicollinearity concerns.
Second, anomalyDetection improves the modeling process to perform multivariate statistical analysis by offering built-in functions to perform Mahalanobis distance, factor analysis, and principal components analysis along with functions to improve the efficiency of extracting and assessing the results from these multivariate approaches.
Third, we demonstrated how the anomalyDetection incorporates the pipe operator (%>%) to allow it to work well in the tidyverse workflow, which helps to improve the overall efficiency of the data analysis process.
It is also important to note that although the authors' focus with this package was to target and improve the analysis of network log-file data, anomalyDetection can also be used for other large data sets that contain arbitrary features that require data aggregation and anomaly analysis. Furthermore, as with any package, we readily admit that further improvements to the package can be made. Future versions of anomalyDetection plan to integrate additional multivariate and time series approaches to offer analysts a wider suite of modeling tools. Integrating plotting functions are also planned for future iterations to further enhance the efficiency of visualizing analytic results. Furthermore, future updates could explore how anomalyDetection could interact with distributable systems such as spark by integrating capabilities from packages such as SparkR. This would improve anomalyDetection's ability to work with Big data architectures.