Introduction

Exploratory data analysis (EDA) is a mission critical task underpinning the predominance of detection development and preparation for cybersecurity-centric machine learning. There are a number of actions that analysts can take to better understand a particular data set and ready it for more robust utilization. In the spirit of toolsmith, and celebration of this being the 150th issue since toolsmith’s inception in late 2006, consider what follows a collection of tools for your security data analytics tool kit.

The University of Maryland’s Center for International and Security Studies (CISSM) Cyber Attacks Database is an ideal candidate for experimental exploration. Per the dataset description, the database “brings together open-source information surrounding a range of publicly acknowledged cyber events on private and public organizations. Events from 2014 through present have been coded to standardize information on threat actor, threat actor country, motive, target, end effects, industry, and country of impact. Source links to the news source are also provided” (Harry & Gallagher, 2018). I asked the project’s principal investigators for data export access as the default UI content is not suited to raw ingestion. The resulting cissm-export.csv (through MAR 2023) file, and all of the code that follows, as well as a Jupyter/Colab notebook for convenient experimentation, are available to you via my CISSM-EDA repository.

We begin with with loading the necessary libraries, data ingestion, data frame construction (tibble), and tsibble creation, a time series tibble. Our series of experiments require correlationfunnel, devtools, forecast, fpp2, CGPfunctions, ggpubr, janitor, tidyverse, tsibble, TTR, and vtree. The following snippet installs packages only if needed. Note that we’re installing dataxray from my fork in order to take advantage of an update I made to the report_xray() function. This update enables the results of a dataxray report to render in a browser automatically, particularly useful when calling the function from a Jupyter/Colab notebook.

my_packages <- c("correlationfunnel", "devtools", "forecast", "fpp2", "CGPfunctions", 
"ggpubr", "janitor", "tidyverse", "tsibble", "TTR", "vtree")                            # Specify your packages
not_installed <- my_packages[!(my_packages %in% installed.packages()[ , "Package"])]    # Extract packages to be installed
if(length(not_installed)) install.packages(not_installed)                               # Install packages
devtools::install_github("holisticinfosec/dataxray")                                    # Install dataxray

With packages installed, we call libraries and build important components for our exercises. Comments are inline for each step.

# Attach the requisite packages

library(dataxray)
library(forecast)
library(fpp2)
library(CGPfunctions)
library(ggpubr)
library(janitor)
library(tidyverse)
library(tsibble)
library(TTR)
library(vtree)

# ingest the CISSM CAD data as a data frame

df <- read_csv("CISSM-export.csv", show_col_types = FALSE)

# shrink the data set to include only event dates and event types

evtType <- tabyl(df, evtDate, event_type)

# convert the reduced data frame to a tibble

df1 <- as_tibble(evtType)

# create an all events tsibble

df1 |>
  mutate(evtDate = yearmonth(evtDate)) |>
  as_tsibble(index = evtDate) -> AllEvents

# create disruptive events tsibble
AllEvents |> select(evtDate,Disruptive) -> disruptive

# create exploitative events tsibble
AllEvents |> select(evtDate,Exploitative) -> exploitative

The disruptive and exploitative variables will be used later as we model, forecast, and plot the results.
We’ll take advantage of the df variable for dataxray, janitor, CGPfunctions, and vtree.

dataxray

dataxray is an interactive table interface for data summaries. It provides an excellent, interactive first look at the data set. I set a specific working directory in my script, set yours as you see fit.

df %>%
  make_xray() %>%
  view_xray()

df %>%
  report_xray(data_name = 'CISSM', study = 'ggplot2')

The result is a columnar view of the CISSM Cyber Attacks Database with variables, observations including missing and distinct data, and interactive figures per variable. The view_xray() function creates an RStudio IDE Viewer pane while report_xray() generates the same result as RMD and HTML files written to your local directory. I’ve hosted the HTML version here if you’d like to interact with it without having to run the analyses yourself.

Figure 1: CISSM CAD dataxray

We learn quickly that the data include 10642 observations and 13 variables.
Of 10642 observations, there are 109 months of data, with 973 unique actors, six actor types, nine motives, four event types and plethora others. Hovering over the interactive figures we discover that motives include espionage, financial, industrial-espionage, personal attack, political-espionage, protest, protest-financial, sabotage, and undetermined. We also learn that event types include disruptive, exploitative, mixed, and undetermined. More on this next where we use janitor to manipulate this data further. dataxray is a great opening salvo in your analysis attack, fulfilling descriptive statistics capabilities admirably.

janitor, CGPfunctions, and vtree

While looking for improved methods to count by group in R I discovered janitor via Sharon’s Infoworld article including “quick and easy ways to count by groups in R, including reports as data frames, graphics, and ggplot graphs” (Machlis, 2020).
janitor includes simple little tools for “examining and cleaning dirty data, including functions that format ugly data.frame column names, isolate duplicate records for further study, and provide quick one- and two-variable tabulations (i.e., frequency tables and crosstabs) that improve on the base R function table().”
The CISSM CAD data sets up perfectly for a variety of table views with counts and/or percentages. While base R’s table() and dplyr’s count() are perfectly useful, a little enrichment never hurt anyone.
As you can see, table and count perform perfectly well.

>table(df$event_type)

  Disruptive Exploitative        Mixed Undetermined 
        3462         5561         1525           94 

>df %>% 
+   count(event_type)
    event_type    n
    
1   Disruptive 3462
2 Exploitative 5561
3        Mixed 1525
4 Undetermined   94

What table and count don’t necessarily provide are the aforementioned enrichment, but tabyl and janitor’s adorn functions fill that void nicely.
First, basic tabyl use yields quick results as seen in the partial output snippet below.

>tabyl(df, event_type, motive)

   event_type Espionage Financial Industrial-Espionage Personal Attack Political-Espionage Protest Protest,Financial Sabotage Undetermined
   Disruptive         0      1335                    0              27                   6     643                 0      240         1211
 Exploitative        25      3077                   86              64                 497     282                 0       42         1488
        Mixed         0      1361                    2               2                  14      41                 1       16           88
 Undetermined         0        57                    0               0                   2       2                 0        1           32

If we chose to adorn our results however the result is all the more beneficial.

>tabyl(df, event_type, motive) %>%
+   adorn_percentages("col") %>%
+   adorn_pct_formatting(digits = 1)

   event_type Espionage Financial Industrial-Espionage Personal Attack Political-Espionage Protest Protest,Financial Sabotage Undetermined
   Disruptive      0.0%     22.9%                 0.0%           29.0%                1.2%   66.4%              0.0%    80.3%        43.0%
 Exploitative    100.0%     52.8%                97.7%           68.8%               95.8%   29.1%              0.0%    14.0%        52.8%
        Mixed      0.0%     23.3%                 2.3%            2.2%                2.7%    4.2%            100.0%     5.4%         3.1%
 Undetermined      0.0%      1.0%                 0.0%            0.0%                0.4%    0.2%              0.0%     0.3%         1.1%

Continuing our journey with event types and motives, a similar construct can be visualized using CGPfunctions.
Per the CGPfunctions description, it’s a package that includes miscellaneous functions useful for teaching statistics as well as actually practicing the art. They typically are not new methods but rather wrappers around either base R or other packages. In this case, we’ll use PlotXTabs2 which wraps around ggplot2 to provide bivariate bar charts for categorical and ordinal data, specifically the event_type and motive variables.

PlotXTabs2(df, event_type, motive, title = "Event Type by Motive")

Figure 2: CISSM CAD PlotXTabs2

PlotXTabs2 stands out as it offers a summary of key frequentist and bayesian information as a subtitle (can be suppressed) as well as a plethora of formatting options courtesy of ggstatsplot (Powell, 2020). Noteworthy in the resulting plot is the fact that across all event types (disruptive, exploitative, mixed, undetermined) the predominant motive was financial. Not surprising, but noteworthy. Regardless, PlotXTabs2 provides an incredibly useful visualization of the CAD dataset with specific variables.

Finally, we conclude this section with use of the vtree package. vtree, or variable trees, displays information about nested subsets of a data frame, in which the subsetting is defined by the values of categorical variables. This is again a perfect option for our subsetted CISSM CAD data. First, we generate a simple vtree based only on event type.

vtree(df, "event_type")

Figure 3: CISSM CAD event type vtree

The result includes a breakdown of event types from the total of 10642 events, including counts and percentage. If you prefer a view without counts you can exclude them as follows; in this case we do so with motives. Note that vtree is incredibly rich with argument features, tailor your visualizations to your liking. Call help("vtree") to learn more.

vtree(df, "motive", showcount = FALSE)

Figure 4: CISSM CAD motive vtree

Finally, we join event types and motives, and pivot our tree to be read vertically.

vtree(df, c("event_type", "motive"), showcount = FALSE, horiz = FALSE)

Figure 5: CISSM CAD vertical event types and motives vtree

The more data you join, the more unwieldy the tree becomes, but nonetheless an intriguing view when zoomed appropriately.

Models

Next, we model the CISSM CAD data for time series forecasting with three well known methods selected for performance with the CISSM CAD dataset. These include naive, SES, and ARIMA.
Naive forecasting uses the most recent observation as the forecast for the next observation.
Simple exponential smoothing (SES) is the method of time series forecasting used with univariate data with no trend and no seasonal pattern.
Autoregressive integrated moving average, or ARIMA, is a statistical analysis model that uses time series data to either better understand the data set or to predict future trends.
Before we initiate the forecasts we use each of the methods to determine which one performs best. For brevity here, we’ll only run the models and plots on exploitative data, but the forecasts_CISSM.R script and the Jupyter/Colab notebook in the GitHub repo run all models and plots with disruptive and exploitative data. Note again that these are subsets of a much broader dataset represented by CISSM CAD. Explore further to your liking.
Important to the modeling that follows, note the evaluation metrics RMSE and MAE. Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) are used to evaluate regression models, tell us how accurate our predictions are, and the amount of deviation from actual values (Acharya, 2021). In essence, the lower the score, the better the performance.

> naive_model_exploitative <- naive(exploitative, h = 12) 
> summary(naive_model_exploitative)

Forecast method: Naive method

Model Information:
Call: naive(y = exploitative, h = 12) 

Residual sd: 20.5334 

Error measures:
                      ME    RMSE      MAE  MPE MAPE      MASE       ACF1
Training set -0.06481481 20.5334 15.41667 -Inf  Inf 0.6228308 -0.3349566

Forecasts:
         Point Forecast      Lo 80     Hi 80      Lo 95     Hi 95
Feb 2023             23  -3.314606  49.31461  -17.24472  63.24472
Mar 2023             23 -14.214473  60.21447  -33.91463  79.91463
Apr 2023             23 -22.578235  68.57824  -46.70590  92.70590
May 2023             23 -29.629213  75.62921  -57.48944 103.48944
Jun 2023             23 -35.841249  81.84125  -66.98992 112.98992
Jul 2023             23 -41.457358  87.45736  -75.57902 121.57902
#snipped

The naive model yields an RMSE score of 20.5 and an MAE score 15.4.

> ses_model_exploitative <- ses(exploitative$Exploitative, h = 12) # RMSE = 18.8, MAE = 13.9
> summary(ses_model_exploitative)

Forecast method: Simple exponential smoothing

Model Information:
Simple exponential smoothing 

Call:
 ses(y = exploitative$Exploitative, h = 12) 

  Smoothing parameters:
    alpha = 0.529 

  Initial states:
    l = 32.3171 

  sigma:  19.0189

     AIC     AICc      BIC 
1157.442 1157.671 1165.517 

Error measures:
                    ME     RMSE      MAE  MPE MAPE      MASE       ACF1
Training set 0.2068036 18.84358 13.92841 -Inf  Inf 0.9034647 0.03522811

Forecasts:
    Point Forecast      Lo 80    Hi 80      Lo 95     Hi 95
110       44.24184 19.8681717 68.61551   6.965531  81.51815
111       44.24184 16.6677635 71.81592   2.070929  86.41275
112       44.24184 13.8020028 74.68168  -2.311874  90.79555
113       44.24184 11.1837445 77.29994  -6.316154  94.79983
114       44.24184  8.7581582 79.72552 -10.025768  98.50945
115       44.24184  6.4880896 81.99559 -13.497539 101.98122
#snipped

The SES model yields an RMSE score of 18.8 and an MAE score of 13.9.

> arima_model_exploitative <- auto.arima(exploitative) 
> summary(arima_model_exploitative)
Series: exploitative 
ARIMA(2,1,1) 

Coefficients:
         ar1     ar2      ma1
      0.4744  0.1671  -0.9653
s.e.  0.1021  0.1007   0.0323

sigma^2 = 343.6:  log likelihood = -467.68
AIC=943.35   AICc=943.74   BIC=954.08

Training set error measures:
                   ME     RMSE     MAE  MPE MAPE      MASE        ACF1
Training set 2.152378 18.19375 13.5603 -Inf  Inf 0.5478338 -0.01933472

Finally, the ARIMA model yields an RMSE score of 18.2 and an MAE score of 13.6. Ultimately, by a small margin, the ARIMA model is most likely to provide the best forecast. Next, we forecast and plot the results. Given that ARIMA is most reliable under these circumstances, we’ll focus on visualizing ARIMA results; you can experiment with naive and SES plots on your own via the scripts or the notebook.

Forecasts & Plots

Generating a plot of AllEvents is as easy as:

autoplot(as.ts(AllEvents))

Figure 6: CISSM CAD AllEvents plot

This is just as easy with disruptive or exploitative events exclusively with the likes of autoplot(as.ts(disruptive)) or autoplot(as.ts(exploitative)).

To forecast the exploitative ARIMA model in an individual plot, utilize:

forecast(arima_model_exploitative) %>% autoplot()

Figure 7: CISSM CAD exploitative events forecast plot

The light and dark areas correspond to the 95% and 80% confidence intervals (CI) respectively.
You can join multiple plots to compare outcomes side by side as follows.

naiveEXP = forecast(naive_model_exploitative) %>% autoplot()
sesEXP = forecast(ses_model_exploitative) %>% autoplot()
arimaEXP = forecast(arima_model_exploitative) %>% autoplot()

multi.pageEXP <- ggarrange(naiveEXP, sesEXP, arimaEXP, nrow = 3, ncol = 1)
multi.pageEXP

Figure 8: CISSM CAD exploitative events multi-model forecast plot

You may be wondering what ARIMA(2,1,1) refers to in our plots. A nonseasonal ARIMA model, which this is, is classified as an “ARIMA(p,d,q)” model, where: p is the number of autoregressive terms, d is the number of nonseasonal differences needed for stationarity, and q is the number of lagged forecast errors in the prediction equation. Therefore, in this case, (2,1,1) is p,d,q found by the auto.arima process indicating that we have two auto-regessive terms, one difference, and one moving average term in our series (Nau, 2020).

Conclusion

Hopefully, this effort has been useful and insightful for security analysts as well as fledgling data scientists in the security realm. It’s no surprise that I orient towards the practices of visualization; I have found all methods deployed here to be useful, effective, and durable for future use. It is my desire that you benefit similarly, and that this opens some doors for you, literally and figuratively.

Cheers…until next time.

References

Acharya, S. (2021, June 15). What are RMSE and Mae? Medium. Retrieved May 3, 2023, from https://towardsdatascience.com/what-are-rmse-and-mae-e405ce230383

Harry, C., & Gallagher, N. (2018). Classifying cyber events. Journal of Information Warfare, 17(3), 17-31.

Machlis, S. (2020, September 10). How to count by group in R. InfoWorld. Retrieved May 1, 2023, from https://www.infoworld.com/article/3573577/how-to-count-by-groups-in-r.html

Nau, R. (2020, August 8). Introduction to ARIMA models. Statistical forecasting: notes on regression and time series analysis. Retrieved May 4, 2023, from https://people.duke.edu/~rnau/411arim.htm#pdq

Powell, C. (2020, November 12). Using PlotXTabs2. The Comprehensive R Archive Network. Retrieved May 2, 2023, from https://cran.r-project.org/web/packages/CGPfunctions/vignettes/Using-PlotXTabs2.html