[1]Department of Computer Science, San Jose State University

The anonymous nature of darknets is commonly exploited for illegal activities. Previous research has employed machine learning and deep learning techniques to automate the detection of darknet traffic in an attempt to block these criminal activities. This research aims to improve darknet traffic detection by assessing Support Vector Machines (SVM), Random Forest (RF), Convolutional Neural Networks (CNN), and Auxiliary-Classifier Generative Adversarial Networks (AC-GAN) for classification of such traffic and the underlying application types. We find that our RF model outperforms the state-of-the-art machine learning techniques used in prior work with the CIC-Darknet2020 dataset. To evaluate the robustness of our RF classifier, we obfuscate select application type classes to simulate realistic adversarial attack scenarios. We demonstrate that our best-performing classifier can be defeated by such attacks,
and we consider ways to deal with such adversarial attacks.

§ INTRODUCTION

Most of us are familiar with the Internet and the World Wide Web (WWW, or web). We regularly access both using web browsers or other networked applications to share information publicly, guided by search engine indexing of the Domain Name System (DNS) over globally bridged Internet Protocol (IP) networks. This publicly accessible and indexed address space is known as the surface web or clearnet. In contrast, the WWW address space which is not indexed by search engines but still publicly accessible is known as the deep web. Private networks within the deep web or networks comprised of unallocated address space are known as darknets and collectively termed the dark web. Figure <ref> illustrates the relationship between these layers of the Internet.

Layers of the Internet <cit.>

The dark web is reached by an overlay network requiring special software, user authorization, or non-standard communication protocols <cit.>. Many darknets afford users anonymity during communication and thus facilitate many criminal activities, including hacking, media piracy, terrorism, trading illegal goods, human trafficking, and child pornography <cit.>.
Researchers are illuminating darknet traffic with machine learning and deep learning techniques, to better identify and inhibit these criminal activities. This research strives to contribute by promoting accurate classification of traffic features from the well-studied CIC-Darknet2020 <cit.> dataset. CIC-Darknet2020 is a collection of traffic features from two darknets, namely The Onion Router (Tor) and a Virtual Private Network (VPN), and also equivalent traffic generated over clearnet sessions using the same applications. We use classic machine learning techniques, such as Support Vector Machines (SVM) and Random Forest (RF) for classification of the network and application type. We also represent traffic features as grayscale images and apply deep learning architectures as classifiers, including Convolutional Neural Networks (CNN) and Auxiliary-Classifier Generative Adversarial Networks (AC-GAN). To assess the issue of extreme class imbalance within CIC-Darknet2020, we explore the option of dataset augmentation, assessing the Synthetic Minority Oversampling Technique (SMOTE) on the performance of our classifiers. Our results show that RF is the most effective at classifying both traffic type and the underlying application types.

Having established baseline classification performance, we consider the confusion of our best classifier in some detail, approaching the problem of darknet traffic detection adversarially. From the perspective of an attacker, we obfuscate the application classes in an attempt to avoid detection. We apply an encoding scheme to transform class features using probability analysis of the CIC-Darknet2020 dataset as a proof-of-concept. We strongly correlate the resulting RF confusion with our obfuscation technique for two attack scenarios, assuming few realistic limitations for traffic modification. We then assess the strength of our obfuscation technique with one defense scenario, by which we demonstrate that we can restore the performance of the RF classifier despite duress. We find that sufficient statistical knowledge of network traffic features can empower either classification or obfuscation tasks.

The rest of this paper is structured as follows. Section <ref> gives a brief background on Tor and VPN and discusses related work on darknet traffic detection. Section <ref> describes the dataset and outlines the methods employed for the experiments. Section <ref> provides background knowledge on the machine learning techniques used in our experiments and gives implementation details. Section <ref> discusses the results of our experiments. Lastly, Section <ref> summarizes our research and explores possible future work.

§ BACKGROUND

In this section, we first discuss the two broad categories of data in our dataset, namely,
Tor and VPN traffic. Then we discuss the most relevant examples of related work.

§.§ The Onion Router

Initially, The Onion Router (Tor) was a project started by the United States Navy to secure government communication. Since 2006, Tor has become a nonprofit with thousands of servers (called relays or relay nodes) run by volunteers across the world <cit.>. Tor clients anonymize their TCP application IP addresses and sessions keys, sending encrypted application traffic through a network of relays <cit.>. An example client application is the Tor Browser, which allows users to browse the web anonymously.
Tor generally selects a relay path of three or more nodes and encrypts the data once for each node using temporary symmetric keys. The encrypted data hops from relay to relay, where each relay node only knows about the previous node and the next node along the path. This design makes it difficult to trace the original identity of Tor clients. Each relay removes a layer of encryption, so that by the last relay, the original data is forwarded to the intended destination as plaintext. Tor then deletes the temporary session keys used for encryption at each node, so that any compromised nodes cannot decrypt old traffic <cit.>.

§.§ Virtual Private Networks

Virtual Private Networks (VPN) are used to ensure communication privacy for individuals or enterprises, and can serve to separate private address spaces from the public Internet. VPN software disguises client IP addresses by tunneling encrypted communications through a trusted server, which acts as a gateway or proxy by routing client traffic to the broader network space. Client data is anonymized behind VPN server credentials before being forwarded to an intended destination, which may be either public or private. Any response traffic is sent back through the VPN server over the encrypted connection for the client to decrypt, ensuring anonymity between the client and recipient. Third parties, such as Internet Service Providers (ISP), will only see the VPN server as the destination of client communications. There are many forms of VPN. Some operate in the network layer, others reside in the transport or application layer <cit.>.

§.§ Related Work

Many researchers have focused on the problem of detecting darknet traffic. However, there are limited public darknet datasets available. The CIC-Darknet2020 dataset used in the experiments reported in this paper was generated by <cit.>. This dataset was also used in prior research, including <cit.>. It has become a well-known darknet traffic dataset due to its accessibility. In their research, <cit.> grouped Tor and VPN together as darknet traffic, while non-Tor and non-VPN were grouped as benign traffic (clearnet). They created $8\times 8$ grayscale images from 61 select features and used CNN to do classification on the dataset. Their CNN model achieved an overall accuracy of 94% classifying traffic as darknet or benign and 86% accuracy classifying the application type used to generate the traffic. The application traffic was broadly labeled as browsing, chat, email, file transfer, P2P, audio streaming, video streaming or VOIP.

Comparatively, <cit.> experimented with classifying traffic and application type by combining a CNN and two other deep-learning techniques: Long Short-Term Memory (LSTM) and Gated Recurrent Units (GRU). They addressed the issue of having an imbalanced dataset by performing Synthetic Minority Oversampling Technique (SMOTE) on Tor, the minority traffic class. They used Principle Component Analysis (PCA), Decision Trees (DT), and Extreme Gradient Boosting (XGB) to extract 20 features before feeding the data into CNN-LSTM and CNN-GRU architectures. Their CNN layer was used to extract features from the input data, while LSTM and GRU did sequence prediction on these features. CNN-LSTM in combination with XGB as the feature selector produced the best F1-scores, achieving 96% classifying traffic type and 89% classifying application type.

The study <cit.> focused on just traffic type from the CIC-Darknet2020 dataset. They used k-Nearest Neighbors (kNN), Multi-layer Perceptron (MLP), RF, DT, and Gradient Boosting (GB) to do binary and multi-class classification. For binary classification, they grouped the data into two classes (benign and darknet), similar to <cit.>. For the multi-class problem, they used the original four classes of traffic type (Tor, non-Tor, VPN or non-VPN). They found that RF was the most proficient classifier for traffic type, yielding F1-scores of 98.7% for binary classification and 98.61% for multi-class classification.

Using the same dataset, <cit.> further broke down the application categories into 11 classes and used Weighted Agnostic Neural Networks (WANN) to classify the data. They proposed automating the tedious design of Artificial Neural Networks (ANN). Unlike regular ANN, WANN do not update neuron weights, but rather update their own network architecture piece-wise. WANN ranks different architectures by performance and complexity, forming new network layers from the highest ranked architecture. Their best WANN model achieved 92.68% accuracy on application layer classification.

The UNB-CIC Tor and non-Tor dataset, also known as ISCXTor2016 <cit.>, was used by <cit.> to classify Tor and non-Tor traffic using Deep Neural Networks (DNN). They built 2 models, DNN-A with 3-layers and DNN-B with 5-layers. DNN-A classified Tor from non-Tor samples with 98.81% accuracy, while DNN-B achieved 99.89% accuracy. For Tor samples, they built a 4-layer DNN to classify 8 application types. This model attained 95.6% accuracy.

In another study, <cit.> generated their own dataset, capturing darknet traffic across 8 application categories (browsing, chat, email, file transfer, P2P, audio, video and VOIP) sourced from 4 different darknets (Tor, I2P, ZeroNet, and Freenet). They used a 3-layer hierarchical approach for classification. The first layer classified traffic as either darknet or normal. In the second layer, samples classified correctly as darknet were then classified by their darknet source. The third layer then classified application type for each of the darknet sources. The techniques <cit.> used for classification include Logistic Regression (LR), RF, MLP, Gradient Boosting Decision Tree (GBDT), Light Gradient Boosting (LightGB), XGB, LSTM, and DT. Their hierarchical method attained 99.42% accuracy in the first layer, 96.85% accuracy in the second layer and 92.46% accuracy in the third layer. Table <ref> provides a summary of the prior work presented in this section.

Summary of previous work
Work     Dataset     Problem considered     Techniques     Results

2*Demertizis, et al. <cit.>
2*CIC-Darknet2020     Only examines     2*WANN     2*92.68% accuracy

11 application types
6*Hu, et al. <cit.>
6*Self-generated     2*Hierarchial approach:     3*LR, RF,     Layer 1:

2*Layer 1: darknet vs clearnet     3*MLP, GBDT,     99.42% accuracy

2*Layer 2: Tor, I2P, ZeroNET       3*LightGB, XGB,     Layer  2:

2*Layer 2:  and FreeNET     3*LSTM, DT     96.85% accuracy

2*Layer 3: 8 application types        Layer 3:

92.46% accuracy
4*Iliadis, et al. <cit.>
4*CIC-Darknet2020     2*Only examines traffic type     3*kNN, MLP     Binary:

2*Binary: darknet vs clearnet     3*RF, DT, GB     98.7% F1-score

2*Multiclass: 4 traffic types        Multiclass:

98.61% F1-score
4*Lashkari, et al. <cit.>
4*CIC-Darknet2020     3*Binary: darknet vs clearnet     4*CNN     Binary:

3*Multiclass: 8 application types        94% accuracy

86% accuracy
4*Sarkar, et al. <cit.>
2*Binary: Tor vs non-Tor     4*DNN      Binary:

2*Multiclass:        99.89% accuracy

2*8 application types within Tor        Multiclass:

95.6% accuracy
4*Sarwar, et al. <cit.>
4*CIC-Darknet2020        3*4 traffic types     3*CNN-LSTM,     Traffic:

3*8 application types     3*CNN-GRU
96% F1-score

89% F1-score

§ METHODOLOGY

A primary goal of this research is to improve upon the state-of-the-art classification of darknet traffic by exploring the performance of SVM, RF, CNN, and AC-GAN as classifiers. We experiment with different levels of SMOTE during a preprocessing phase, oversampling the minority classes of the CIC-Darknet2020 dataset to assess the effects of data augmentation and class balance on classifier performance. We also experiment with representations of the darknet traffic features as 2-dimensional grayscale images for CNN and AC-GAN. Then we test the robustness of our best-performing classifier through a few obfuscation scenarios which serve to simulate adversarial attacks, assuming both the perspectives of an attacker and defender.

We apply statistical knowledge of the dataset to obfuscate our data features, disguising one or more classes as others. We explore three scenarios whereby we either obfuscate the training data, the test data or both. Obfuscating just the test data simulates an attack scenario in which traffic data is disguised while our classifier is yet unaware of the attack, and thus we can only apply previously trained models without a chance to learn of the obfuscation. Obfuscating just the training data simulates a scenario in which an attacker has accessed our training data to poison it, such that we train our classifier with malformed assumptions or outright malicious supervision. A third scenario supposes we collect some of the obfuscated traffic data before training our classifier, and thus have a chance to update our classification models to detect obfuscated test data.

§.§ Dataset

The CIC-Darknet2020 dataset <cit.> is an amalgamation of two public datasets from the University of New Brunswick. It combines the ISCXTor2016 and ISCXVPN2016 datasets, which capture real-time traffic using Wireshark and TCPdump <cit.>. CICFlowMeter <cit.> is used to generate CIC-Darknet2020 dataset features from these traffic samples. Each CIC-Darknet2020 sample consists of traffic features extracted in this manner from raw traffic packet capture sessions. CIC-Darknet2020 consists of 158,659 hierarchically labeled samples. The top level traffic category labels consist of Tor, non-Tor, VPN or non-VPN. Within those categories, samples are further categorized by the types of applications used to generate the traffic. These subcategories include audio-streaming, browsing, chat, email, file transfer, P2P, video-streaming and VOIP. Table <ref> details the applications that are used to generate each type of traffic at the application level.

CIC-Darknet2020 application classes <cit.>
Application class     Applications considered
Audio-Streaming     Vimeo and YouTube

Browsing     Firefox and Chrome

Chat     ICQ, AIM, Skype, Facebook and Hangouts

Email     SMTPS, POP3S and IMAPS

File Transfer     Skype and FileZilla

P2P     uTorrent and Transmission (BitTorrent)

Video-Streaming     Vimeo and YouTube

VOIP     Facebook, Skype and Hangouts

§.§ Preprocessing

The CIC-Darknet2020 dataset has samples with missing data, more specifically feature values of . We remove samples with these values to improve the performance of the classifiers. As shown in Table <ref>, there are significantly less Tor samples compared to the other traffic categories. Prior work using this dataset eliminated CICFlowMeter flow labels, including: Flow Id, Timestamp, Source IP and Destination IP. The Flow Id and Timestamp are eliminated for this research as well. However, to retain as much information as possible from the CIC-Darknet2020 dataset, we separate each octet of the source and destination IP addresses into their own feature columns. Preliminary tests run on the dataset with and without these IP octet features indicate an improvement in the performance of the classifiers when the

IP information is retained.
Thus our feature dataset contains 72 features total after this preprocessing step.

Samples per traffic category
Traffic Category     Samples
Non-Tor     93,357

Non-VPN     23,864

Tor     1,393

VPN     22,920

The CIC-Darknet2020 dataset was scaled by min-max normalization, which applies the equation

\mbox{{\tt normalizedValue}} = \frac{(\mbox{{\tt value}} - \mbox{{\tt min}})}{(\mbox{{\tt max}} - \mbox{{\tt min}})}

to every value in each feature column to scale the feature values between 0 and 1. We also apply min-max normalization to our new IP octet feature columns.

§.§.§ Data Balancing

The CIC-Darknet2020 dataset does not have balanced sample counts among traffic and application classes, as shown in Tables <ref> and <ref>. To explore the effect of reducing this imbalance on the classification task, we oversample each minority class using SMOTE. SMOTE interpolates linearly between feature values to produce new samples <cit.>. We experiment with different levels of oversampling
including 0% (no SMOTE), 20%, 40%, 60%, 80% (partial SMOTE), and 100% (full SMOTE). SMOTE is performed on all classes with less than the oversampling threshold as compared to the class with the largest sample count. Note that 100% SMOTE results in an equal number of samples for each class, while lower thresholds of SMOTE result in an equal number of samples among only those classes which are oversampled.

Samples per application category
Application Category     Samples
Audio-Streaming     18,065

Browsing     32,809

Chat     11,479

Email     6,146

File Transfer     11,183

P2P     48,521

Video-Streaming     9,768

VOIP     3,567

§.§.§ Data Representation

SVM and RF both use each dataset sample in its original format, which is a 1-dimensional array. However, we reshape each sample to be 2-dimensional for CNN and AC-GAN. Intuitively, the data is reshaped as $9\times9$ grayscale images, where each of our 72 features is represented as a single pixel with the remaining pixels produced by zero padding. The pixels are ordered as their respective features appeared in the CIC-Darknet2020 dataset, starting at the top left corner of the image as shown in Figure <ref>,

where each row represents samples from an application class, color-coded for readability.

Data as 2-D images in original order

Both CNN and AC-GAN convolve local structures within the 2-D images, so adjacent pixels
play an important role in classification. Therefore, we experiment with strategies to reorder the data to achieve better performance. We order the pixels by RF feature importance, starting at the top left corner of the image, and also reorganize the pixels spiraling outward from the center of the image. This latter strategy tends to group pixels with larger values toward the center of each image, as shown in Figure <ref>.

Data in 2-D sorted by RF feature importance and centered

§.§.§ Data Augmentation Experiment

We experimented with AC-GAN as an alternative to SMOTE, to generate realistic artificial samples in an attempt to augment our dataset as another way to address the issue of class imbalance. However, we abandoned this approach as we found that the fake images generated by AC-GAN are consistently detectable by a CNN model with near perfect accuracy,
ranging from 99% to 100% accuracy within a few epochs of training. The depth of our neural network architecture for AC-GAN was constrained by the input image size of our data. Therefore, we were unsuccessful in our attempt to use AC-GAN to augment our data. An example of some fake samples compared to real samples can be found in Figure <ref>.

Chat class (4 fake and 4 real examples)

§.§ Evaluation Metrics

In our experiments, we use
accuracy and F1-score to measure the performance of each classifier. Accuracy is computed as the total number of correct predictions over the number of samples tested. The F1-score is the weighted average of precision and recall metrics, which is better for unbalanced datasets like CIC-Darknet2020. Similar to accuracy, F1-scores fall between 0 and 1, with 1 being the best possible score. F1-score is computed as

\mbox{{\tt F1}} = {\tt 2} \times \frac{(\mbox{{\tt Precision}} \times \mbox{{\tt Recall}})}{(\mbox{{\tt Precision}} + \mbox{{\tt Recall}})}

Precision calculates the ratio of samples classified correctly for the positive class, while recall measures the total number of positive samples that were classified correctly. Precision and recall are computed as

\mbox{{\tt Precision}} = \frac{\mbox{{\tt True Positives}}}{(\mbox{{\tt True Positives}} + \mbox{{\tt False Positives}})}

\mbox{{\tt Recall}} = \frac{\mbox{{\tt True Positives}}}{(\mbox{{\tt True Positives}} + \mbox{{\tt False Negatives}})}

§ IMPLEMENTATION

This section details the implementation of the experiments that
we mentioned in Section <ref>.
All experiments are coded in Python. The Imblearn library is used to implement SMOTE to balance the dataset, while the package Scikit-learn is employed to run SVM and RF experiments. The Tensorflow and Keras libraries are utilized to implement CNN and AC-GAN. From the Scikit-learn library, the metrics module is used to evaluate the F1-scores and accuracy of the classifiers and the StratifiedKFold function is applied to perform 5-fold cross validation. Graphs are generated with the Matplotlib and Seaborn libraries.

All experiments in this research are executed on one of two personal computers, as detailed in
Table <ref>. We exploit a graphics processing unit (GPU) in the second computer to decrease the runtime duration of our more computationally demanding experiments, that is those using neural nets to process 2-D image representations. However, for portability and convenience we also use the first computer, which is a light-weight laptop.

Computer hardware
Computer Processors     Experiments

2*CPU: 8-core Intel(R) CORE(TM)        SVM

2*i7-8550U @ 1.80GHz              RF

CPU: 12-core Intel(R) Xeon(R)        2*CNN

W-10855M @ 2.80GHz              2*AC-GAN

GPU: NVIDIA Quadro           2*SMOTE

RTX 5000

§.§ Background of Classification Techniques

This section describes the machine learning and deep learning concepts that we
apply to classification in our experiments. These include SVM, RF, CNN, and AC-GAN.

§.§.§ Support Vector Machines

SVM are supervised machine learning models frequently used for classification. SVM find one or more hyperplanes to separate labeled training data while maximizing the margin of decision boundaries between classes. Maximizing the decision margin minimizes the classification error. The data must be vectorized into linear feature sets, but non-linear data can also be encoded with some success. Scaling the feature values across training samples allows coefficients of the hyperplanes (weights) to be ranked by relative importance. SVM rely on the kernel trick to map data into a higher dimensional space. The idea behind the kernel trick is that a larger feature space makes it easier to find hyperplanes between classes. Once hyperplanes are found, SVM can classify new samples into one of the classes as separated by the hyperplanes <cit.>. For our research, we perform preliminary tests to determine the best kernel for our dataset, with the result being the Gaussian radial basis function (RBF).

§.§.§ Random Forest

RF is an ensemble method that generalizes DT. While DT is a simple and efficient classification algorithm, it is sensitive to variance in the training data and prone to overfitting. RF compensates for these deficiencies by generating many subsets of the dataset, then randomly selecting features (with replacement) to train a DT for each subset. This process is called bootstrapping. To classify, RF takes the majority vote from all resulting DT in a process called aggregation. Together bootstrapping and aggregation is referred to as bagging <cit.>. RF also produces a useful attribute called feature importance, which is a list of all the features ranked by the mean entropy within the DT. Feature importance tells us how influential each feature is when classifying samples with RF.

§.§.§ Convolutional Neural Networks

CNN are a unique type of neural network that operate on local structures, making them ideal for image analysis. CNN are composed of an image input layer, convolution and pooling layers and an output layer that produces a single vector of class scores. The neural nodes of convolution layers are connected to only a small portion of the nodes of the previous layer, with the exception of dense layers which have connections to every node. Convolution, pooling and dense layers are the fundamental components of any CNN architecture. In convolution layers, an input image is convolved with initially randomized filters to produce local structure maps that are joined to create the output of this layer. The filter windows slide across the input image, computing a dot product of the filter and image at each location. Pooling layers decrease total training time by reducing the dimensionality of the resulting feature maps, concentrating effort on the highest-value features <cit.>. There are different types of pooling; for this research, we use max pooling.

Our CNN architecture is based on that described in <cit.>. We experiment with hyperparameters, running different combinations where we adjust the following.

* Initial number of convolution filters (9, 32, 64, 81)

* Filter size ($2\times2$, $3\times3$)

* Percentage dropout $(0.2, 0.5)$

* Number of nodes in the first dense layer (72, 256)

All these architectures yield accuracies within the range of 86% to 88% when classifying application type. Therefore, we select the architecture that produces the highest accuracy.

We use Adam for our optimizer and sparse categorical cross entropy for our loss function.

Dropout is a common technique used in neural networks with fully-connected layers to combat overfitting. However, it is found to be not as effective with convolution layers. A better suited regularization technique for CNN is to cut out sections of the input images. Cutouts delete features to force CNN to learn from the other parts of each image during training. It is comparative in effect to dropout except that it operates on the input stage rather than the intermediate layers <cit.>. We implement cutouts by creating feature masks of equivalent size to our input image. We experiment with different cutout sizes including $2\times2$, $3\times3$, and $4\times4$ and randomize the position of the cutout within the mask. The masks are applied to our training data with matrix multiplication. Refer to Figure <ref> for some examples of masks with $3\times3$ cutouts.

Examples $9\times9$ images with $3\times3$ cutouts

§.§.§ Auxiliary-Classifier Generative Adversarial Network

Generative Adversarial Networks (GAN) are comprised of two neural network architectures, a generator and a discriminator, that compete in a zero-sum game during training. The generator takes noise from a latent space as input and produces images that feed into the discriminator. The discriminator is given both real and generated images and is tasked to classify them as either real or fake. The discriminator error is then fed back into the generator as direction for image generation. AC-GAN is an extension of this base GAN architecture, taking a class label as additional input to the generator while predicting this label as part of the discriminator output. The objective of the AC-GAN generator is to minimize the ability of the discriminator to distinguish between real and fake images but also maximize the accuracy of the discriminator when predicting the class label <cit.>. Besides using the AC-GAN generator in our data augmentation experiment, we also explore the secondary class prediction output of the discriminator as a classifier.

Our AC-GAN architecture is inspired by the ImageNet model described in <cit.>. However, since that architecture was built for image sizes $32\times32$ and bigger, we modify that architecture to accommodate our $9\times9$ image size by reducing the number of convolution and transposed convolution layers in the discriminator and generator, respectively.

We fine-tune our hyperparameters by experimenting with the following.

* Latent space size (81, 100)

* Initial number of convolution filters (15, 40, 64, 192, 202, 384, 500, 1500)

* Number of nodes in the first dense layer (31, 81, 128, 384, 405, 768, 1000, 3000)

* Filter size ($3\times3$, $5\times5$)

* Stride size ($2\times2$, $3\times3$)

We observe accuracies within the range of 70% to 73% when classifying application type with these hyperparameters. The best-performing architecture with the shortest runtime duration is used in this research. Tables <ref> and <ref> detail our generator and discriminator architecture, respectively.

We feed training data to our AC-GAN model in batches of 64 samples. Batch normalization (BN) layers are applied between convolution layers to regularize the training gradient step size. This is thought to smooth local optimization steps and stabilize training, thereby accelerating convergence of GAN models <cit.>.

AC-GAN generator architecture
Layer Operation     Kernel     Strides     Depth     BN     Activation

$1\times 1\times 100$ Input A (Latent Space)

Dense A                 405           ReLU

$1\times 1\times 1$ Input B (Feature Noise)

$8\times 32$ Class Embedding for B                 256

Dense B                 1

Merge A + B                 406

Conv2DTranspose     $5\times 5$     $3\times 3$     202         ReLU

Conv2DTranspose     $5\times 5$     $3\times 3$     1           Tanh

AC-GAN discriminator architecture
Layer Operation     Kernel     Strides     Depth     BN     Dropout     Activation

9x9x1 Input (Image)

Conv2D     $3\times3$     $2\times2$     32          0.5     Leaky ReLU

Conv2D     $3\times3$     $1\times1$     64         0.5     Leaky ReLU

Conv2D     $3\times3$     $2\times2$     128         0.5     Leaky ReLU

Conv2D     $3\times3$     $1\times1$     256         0.5     Leaky ReLU

Dense             1             Sigmoid

Dense             8             Softmax
Leaky ReLU Slope                         0.2

Weight Initialization                         Gaussian $(\sigma=0.02)$

Optimizer                         Adam $(\alpha=0.0002,  \beta_1=0.5)$

§.§ Adversarial Attacks

Our adversarial attacks rely on on obfuscation, which serves to disguise application classes
based on applied probability analysis. We select
application classes to disguise as other classes based on minimum and maximum sum statistical distance between all class features, as specified in Algorithm <ref>.

Class feature probability distributions
[1]
bins = some discrete bins partitioning values 0 to 1
We use 100 bins.
A = classA feature probability distributions
B = classB feature probability distributions
classDistance = 0
each distributionA, distributionB in A, B
featureDistance = cdist(distributionA, distributionB)
classDistance += featureDistance
Manhattan Sum

We also select a third class transformation to perform based on maximal classifier confusion, whose sum statistical distance between class features is notably low but not the minimum between classes. We ensure our class transformation can be decoded by encoding features with a deterministic algorithm, Algorithm <ref>, but impose no other restrictions on feature transformation.

Disguise one class sample as another class sample
[1]
To decode, reverse A and B
bins = some discrete bins partitioning values 0 to 1
We use 100 bins.
A = classA feature probability distributions
B = classB feature probability distributions
each featureValue at featureIndex in the sample
featureBin = the bin which contains featureValue
DCPD = A[featureIndex] - B[featureIndex]
AtoB = sorted DCPD from maximum to minimum
BtoA = sorted DCPD from minimum to maximum
oldBin = where AtoB[oldBin] == featureValue
Red arrows in Figure <ref>
newBin = BtoA[oldBin]
Black arrows in Figure <ref>
newValue = featureValue - (bins[oldBin] - bins[newBin])
sample[featureIndex] = newValue

We start by generating normalized histograms of feature values per class to assess the probability at which values occur within each class over our entire dataset. To decide which classes to obfuscate, we examine the sums of the distances between feature probability distributions from each class to each other class, depicted as an $8\times8$ table in the case of application type as provided in Table <ref>. We use the cdist function of the scipy Python library to calculate the Euclidean distance between probability distributions. This allows us an estimation of the overall difference between classes while considering all feature probability distributions equally. We pick the classes with the minimum and maximum sum of statistical distances between features, changing class 0 (audiostreaming) to class 5 (P2P) and class 3 (email) to class 7 (VOIP). We also examine the confusion matrix for our best-performing classifier—RF, shown in Figure <ref>. RF is observed to be most confused between class 2 (chat) and class 3 (email), so we decide to additionally obfuscate class 2 with class 3. We arbitrarily choose to transform lower numbered classes to higher numbered classes (e.g. disguising class 2 as
class 3 instead of class 3 as class 2).

Statistical distances between pairs of application classes

Class     0     1     2     3     4     5     6     7

0     0    	25.129     18.709     21.041     23.656    	28.195 [gray].9     18.371     21.903
1     25.129     0     23.518     21.958     12.884     12.098     16.728     23.623
2     18.709     23.518     0     9.841 [gray].75     22.613     25.294     18.408     9.901
3     21.041     21.958     9.841 [gray].75     0     21.51     23.021     18.031     6.859 [gray].6
4     23.656     12.884	    22.613     21.51     0     15.651     14.605     23.211
5     28.195 [gray].9     12.098     25.294     23.021     15.651     0     21.085     24.451
6     18.371     16.728     18.408     18.031     14.605     21.085     0     20.089
7     21.903     23.623     9.901     6.859 [gray].6     23.211     24.451     20.089     0

Confusion matrix for best RF classifying level 2

Our obfuscation algorithm first calculates the difference in class probability distributions (DCPD) for each feature between 2 classes, denoted as A and B, and sorts each distribution from maximum to minimum. Intuitively the index of each DCPD maximum corresponds to each feature value most probably belonging to the positive class A while the minimum corresponds to each feature value most probably belonging to the negative class B. To obfuscate a sample, we then transform individual feature values by subtracting the difference in bin thresholds between the original feature bin and a target bin for obfuscation. To choose target bins for this transformation, we create a 1-to-1 map of the sorted indices of each DCPD with a reverse sort of the same DCPD. This ensures a transformed sample feature could be decoded later given the feature DCPD for a class obfuscation vector.

An example visualization of DCPD bin mapping for the transformation of the most
common feature 0 values from class 2 to class 3 is provided in Figure <ref>. Reversing the 1:1 bin map facilitates decoding of obfuscated class sample feature values back to their original values. To do this we add back the same difference in bin thresholds which we subtracted earlier, thus applying each feature DCPD between known classes as a decoder key to undo an expected class obfuscation for a particular feature. To test this method of class obfuscation, we performed the three experiments described in Section <ref> and summarized in Table <ref>, with RF as the classifier.

Obfuscation scenarios
1r|     2c|What is obfuscated?      2*Scenario description

Training data     Test data
3*Case 1                  3*       Simulates a novel attack

3*                        3*       where we apply an outdated

3*                        3*       model for classification
3*Case 2               3*        Simulates an attack

3*                        3*       on our training data,

3*                        3*       poisoning the classifier
3*Case 3               3*     3*       Simulates a novel defense

3*                        3*       where we train our model

3*                        3*       on some obfuscated data

§.§.§ An Obfuscation Example

To better illustrate Algorithm <ref>, we will walk through a simple example where we are given a sample from class 2 and we want to transform this sample to look like class 3. Let us start with the first feature, feature 0. We note the value of this feature for class 2; call this value $v$.

Suppose, for example, that $v = 0.178142$. We allocate 100 equal-width bins ranging from 0 to 1, so that bin $b_0$ corresponds to values 0.00 to 0.01 and so on. Given the value of $v$, we find the bin that $v$ falls into. The value $v= 0.178142$ is in bin $b_{17}$, which contains values between 0.17 to 0.18. Bin $b_{17}$ is indicated by the red arrows in Figure <ref>. We then flip the sorted DCPD index at $b_{17}$ to locate our target bin, indicated by the black arrows in Figure <ref>. This target bin $b_{58}$, which contains values between 0.58 to 0.59. To obfuscate, we subtract the difference between $b_{17}$ and $b_{58}$ from $v$. In this example, our new transformed value is

\widetilde{v} = 0.178142 - (0.17 - 0.58) = 0.588142

which falls into the target bin $b_{58}$. We repeat this for all the features to transform the sample from class 2 to class 3.

Visualization of obfuscation example

Note that this obfuscation technique is designed to maximize the effectiveness of
a simulated adversarial attack. Our approach ignores practical limitations on the ability
of attackers to modify the statistics of the data. Hence these simulated attacks
can be considered worst-case scenarios, from the perspective of
detecting darknet traffic under adversarial attack.

§ RESULTS AND DISCUSSION

In this section, we consider a wide range of experiments. First, we determine
which of the three 2-D image representation techniques discussed in
Section <ref> is most effective. Then we consider the use
of cutouts, which can serve to reduce overfitting and improve accuracy in CNNs.
We then turn our attention to the imbalance problem, with a series of SMOTE
experiments. We conclude this section with an extensive set of experiments
involving various adversarial attack scenarios.

§.§ Data Representation Experiments

We evaluate CNN and the AC-GAN discriminator given different 2-D pixel representations of the data features. All of our 2-D representations of the data are of size $9\times9$, where each pixel is a feature. The pixels in the original representation follow the order that the features appear in the CIC-Darknet2020 dataset. We hypothesize that grouping the pixels together would have a positive effect on the performance of our classifiers since convolutions operate on local structures. Our results show that CNN performs best when the pixels are sorted by RF feature importance and then grouped together at center of the image. However, this is not true for the AC-GAN discriminator. AC-GAN does better using the original data representation, contrary to our hypothesis. Table <ref> shows the results for these experiments.

2-D data representation performance for application classification
1r|     2c|CNN     2cAC-GAN

Accuracy     F1-scores     Accuracy     F1-score

Original     0.889     0.887     0.753     0.738

Shaped with RF feature importance     0.890     0.887     0.753     0.731

Shaped with RF feature importance & centered     0.891     0.889     0.742     0.729

§.§ Cutout Experiments

Initially, our CNN model is able to achieve 88% accuracy classifying application type within 15 epochs. However, we notice that overfitting starts to occur the longer we run our model. To reduce overfitting, we apply cutouts to the training data. We experiment with different cutout sizes: $2\times2$, $3\times3$, and $4\times4$. We observe that cutouts allow our CNN to train for a longer period of time without overfitting. The loss graphs in Figure <ref> show how the CNN model overfits after 20 epochs in the original execution but does not overfit with cutouts. There is little difference in the effects of applying $2\times2$ compared to $3\times3$ cutouts. Both delay overfitting at the same rate and the accuracies for both linger at 88%. Notably, we witness a 1% decrease in accuracy with $4\times4$ cutouts. As our images are only $9\times9$ pixels, a $4\times4$ cutout may delete too much information from the image, negatively affecting the accuracy. While cutouts address the issue of overfitting, we find that more training does not significantly improve the performance of CNN for our dataset.

The effects of cutouts on overfitting for CNN

§.§ SMOTE Experiments

We compare the performance of our classifiers with various levels of SMOTE, performing SMOTE to oversample the training data before training each classifier to classify traffic and application types; the results from these experiments appear in Tables <ref>
and <ref>, respectively.
We observe
that reducing class imbalance this way does not have a significant effect on the performance of RF and CNN. In contrast, we do see improved performance with SVM for both traffic and application classification tasks as the level of applied SMOTE increases. We also notice improved performance with increasing SMOTE level for AC-GAN when doing traffic type classification. We observe that partial SMOTE improves the performance of AC-GAN for application type classification, while noting that the performance decreases with full SMOTE.

Traffic classification F1-scores for increasing levels of SMOTE
SMOTE level     SVM     RF     CNN     AC-GAN

0%    0.986      0.998      0.998     0.974

20%    0.993      0.998      0.995     0.980

40%    0.993      0.998      0.995     0.984

60%    0.993      0.998      0.995     0.986

80%    0.993      0.998      0.996     0.987

100%    0.993      0.998      0.995     0.987

Application classification F1-scores for increasing levels of SMOTE
SMOTE level     SVM     RF     CNN     AC-GAN

0%    0.834      0.922      0.887     0.738

20%    0.839      0.920      0.883     0.750

40%    0.842      0.921      0.883     0.762

60%    0.846      0.921      0.887     0.768

80%    0.847      0.920      0.888     0.767

100%    0.848      0.920      0.885     0.759

Our RF model without SMOTE outperforms the state-of-the-art F1-scores for both traffic and application classification tasks. We observe a 1.1% improvement for traffic classification as compared to <cit.>, where they also found RF to be their best classifier. Study <cit.> only classified traffic type, thus no application type performance is available for comparison. For application classification, our RF model achieved a 3.2% increase over <cit.>. In addition, our CNN model outperformed
the CNN results in <cit.> by 2.8%. We are only able to compare classification results for application type with study <cit.> because they approach traffic type classification as a binary problem while we address it as a multiclass problem. Table <ref> summarizes the best performance of our classifiers with relevant prior work. Overall, RF is our best-performing classifier and AC-GAN discriminator was shown to be our worst-performing classifier at any SMOTE level.

Best F1-score compared to prior work
2*Data type
3c|Previous work     4cOur results

CNN-LSTM <cit.>     RF <cit.>     CNN <cit.>     SVM     RF     CNN     AC-GAN

Traffic     0.960     0.987     —     0.993     0.998     0.998     0.987

Application     0.890     —     0.860     0.848     0.922     0.888     0.759

§.§ Adversarial Attack Experiments

With continual improvement in the accuracy of darknet traffic detection by machine learning and deep learning techniques, it is realistic to anticipate that attackers will attempt to find ways to circumvent detection by modifying the profile of their application traffic. For example, someone pirating copyrighted media with P2P applications might disguise their illegal activity as VOIP traffic to avoid prosecution. We show obfuscation of traffic in this fashion can be accomplished by modifying traffic feature values, understanding that this process is most feasible and desirable at the application layer. Also, if an attacker were to discover the methods we use for classification and pollute our training data, then our classifiers could be compromised, allowing the attacker to avoid detection without modifying any of their traffic features.

For this experiment we assume the role of an attacker on the network, with the goal of modifying traffic features such that classes are incorrectly classified or entirely undetected. This could represent covert illegal activity that an attacker wishes to hinder the detection of, common examples being P2P or file-transfer applications. Realistically, traffic features common to one application class could be modified at the application layer to appear more similar to other application classes. An attacker could do this by writing a custom overlay application to change various features, such as the number of packets sent, their communication intervals, port assignment or many others. In our experiments, we disguise class 0 as class 5 (originally the most different), class 2 as class 3 (the classes which most confused our RF classifier)
and class 3 as class 7 (originally the most similar).

In scenario case 1, we train our RF classifier on the original application class data, then test the same model with an obfuscated class in the test dataset. This represents a hypothetical scenario where an attacker modifies the traffic features of one class at the application layer, perhaps with an overlay application. We demonstrate that our method of obfuscation is able to defeat our best classifier in this scenario, significantly reducing detection of the obfuscated class, as well as overall classifier accuracy. Before obfuscation, RF classifies application classes with an accuracy of 92.3% without SMOTE. After obfuscation of our three class choices mentioned in the previous paragraph, the overall RF accuracy decreases to 80.8%, 85.4%, and 88.7% respectively, for application classification without SMOTE.

The confusion matrices in Figure <ref> show that RF consistently misclassifies each class we obfuscate, actually detecting no samples in the case of an obfuscated class 3, where the dashed circle indicates the class we are obfuscating and the solid circle indicates the class we intended it to appear as.
However, RF did not misclassify classes 0 and 3 as the expected classes 5 and 7, respectively. Instead, the confusion matrices (a) and (c) in Figure <ref> reveal that RF mostly categorizes class 0 and 3 as class 6 and 2, respectively. It may be relevant that our obfuscation method does not account for any interdependance between traffic feature values, obfuscating each feature independently.

(a) class $0 \longrightarrow 5$

(b) class $2 \longrightarrow 3$

2c(c) class $3 \longrightarrow 7$

Confusion matrices for obfuscation scenario case 1

In scenario case 2, we train our RF classifier with an obfuscated class in the training dataset, then test the model with the original application class data. We consider a hypothetical scenarios where an attacker entirely poisons our training data, perhaps by injecting malware into our database or by intercepting our traffic capture data stream. We find the attacker could prevent an entire class from being predicted by our best classifier when the training data for a class is entirely obfuscated. We see this trend in the all three confusion matrices in Figure <ref>,
where in each case, 	the dashed circle indicates the class we are obfuscating and the solid circle indicates the class we intended it to appear as.
Notice that the entire row in the confusion matrix is zeroed out, indicating that the class was never predicted for classification by RF. Similar to scenario case 1, the overall RF accuracy
decreases to 82.0%, 86.5%, and 89.4% respectively, for application classification without SMOTE. The decrease in overall accuracy for each class obfuscation again tracks with the count of samples within application classes 0, 2, and 3. However, as the obfuscated class is never considered for prediction by RF, we observe a lesser overall accuracy decrease as compared to scenario case 1.

(a) class $0 \longrightarrow 5$

(b) class $2 \longrightarrow 3$

2c(c) class $3 \longrightarrow 7$

Confusion matrices for obfuscation scenario case 2

In scenario case 3, we train our RF classifier with the same obfuscated class in both the training dataset and the test dataset. We obfuscate only a small portion of the training data while still obfuscating all of the test data for each of class 0, 2, and 3. We experiment with the percentage of training data we obfuscate. This represents a hypothetical scenario where the obfuscation algorithm has been obfuscating network traffic long enough to pollute a small portion of a network traffic population. A defender then updates the classifier to include this small portion of obfuscated class data at training time, with increasing exposure to the obfuscated data over time. As our dataset is split into 80% training data and 20% test data, we decide to limit the training dataset exposure of obfuscated class data to 20% of the total training dataset. We choose to decrement this value logarithmically with three total sub-scenarios representing 0.2%, 2%, and 20% obfuscation exposure, expecting that with more exposure to the obfuscated class data, our classifier will adapt and outperform the obfuscation algorithm to correctly classify the obfuscated class in our test dataset.

We find that 20% exposure of our obfuscation algorithm to the RF training data is sufficient for RF to predict the disguised classes with high accuracy, defeating our obfuscation technique as shown in Table <ref>. Note that the overall accuracies reported for scenario case 3 are higher than our RF benchmark score of 92.2%. We modify the test dataset in both scenario cases 1
and 3, so the resulting accuracies of those scenarios cannot be compared to the results of prior work. Our results with lower exposure levels of 2% and 0.2% reveal a trend. Class 0 appears to be the most difficult class for our algorithm to obfuscate, while class 3 appears to be the easiest for our algorithm to obfuscate. Class 2 is somewhere in between, providing a loose correlation to our metric of statistical distance between classes and the performance of our obfuscation algorithm. We observe this trend in Table <ref> under “Class Accuracy” for “Case 3.”

Class accuracies for obfuscation scenarios
2*Case     2*Obfuscation      3c||Overall Accuracy     3cClass Accuracy

(0,5)     (2,3)     (3,7)     0     2     3

Original     —     —     —     —     0.907     0.846     0.821

1     —     0.808     0.854     0.887     0.100     0.006     0.000

2     —     0.820     0.865     0.894     0.000     0.000     0.000

3     20.0%     0.947     0.946     0.939     0.998     0.993     0.997

3     2.0%     0.935     0.939     0.891     0.958     0.921     0.120

3     0.2%     0.820     0.859     0.887     0.503     0.247     0.000

§ CONCLUSION AND FUTURE WORK

In this research, we classified the CIC-Darknet2020 network traffic samples and assessed SVM, RF, CNN, and
the discriminator model of AC-GAN as our classifiers. We classified 4 traffic and 8 application classes while fine tuning the classifier hyperparameters. We experimented with different levels of SMOTE to assess class imbalance in the dataset and explored 2-D representations of the traffic features for CNN and AC-GAN. We also approached the issue of darknet detection adversarially, from the perspective of an attacker hoping to confuse our best classifier. We demonstrated that we could effectively obfuscate application class traffic features. We then correlated the underlying statistics of the CIC-Darknet2020 dataset to the performance of this algorithm assuming some hypothetical scenarios for added realism.

Among the tested machine learning classifiers, RF was found to be the most proficient at classifying darknet traffic for both traffic and application types. It yielded 99.9% F1-score for traffic classification and 92.3% F1-score for application classification, outperforming the state-of-the-art studies on CIC-Darknet2020 <cit.>. Figure <ref>
provides a visual comparison of our best results with those of prior work.

Classification accuracy compared to previous work

Our research was limited by the availability of darknet traffic datasets. We selected the CIC-Darknet2020 dataset because it is frequently cited and publicly accessible, however the dataset suffers from a substantial imbalance. We attempted to compensate for this class imbalance by generating artificial samples with AC-GAN and SMOTE. The artificial SMOTE samples did not appear realistic for this dataset, which contains primarily non-linear feature values. Seeking to improve the quality of artificial samples, we assessed AC-GAN as a sample generator. However, our AC-GAN generated samples were not useful for data augmentation purposes. An approach future research might consider is to use clustering to group samples within a class, then train one GAN per cluster to generate samples. Other variations of GAN might also be better suited for multiclass sample generation and could generate more authentic samples.

We kept our obfuscation fairly basic, with the goal being to simply demonstrate that we could confuse our best classifier, with few restrictions imposed on the algorithm. Under realistic scenarios, it may be impossible to so easily modify features which define darknets such as Tor and VPN, but possible to obfuscate traffic features at the application layer such as those produced by CICFlowMeter analysis. We introduced a loose correlation to one statistical metric, an independent (Manhattan) sum of distances between DCPD across all sample features. We noted that 2 out of the 3 classes we chose to obfuscate were misclassified not as the intended classes, but with a majority of predictions distributed among other classes. This results from the fact that our obfuscation metric does not account for the statistical relationship between more than two classes, nor does it account for any dependency between the CIC-Darknet2020 feature values.

There is much remaining that could be done to extend the adversarial obfuscation analysis presented this paper. Real traffic features could be modified on live network traffic (e.g. changing IP addresses, ports, packet lengths or intervals), or select features could be prohibited from modification during obfuscation, which is likely to be a realistic constraint. An even larger task is to explore the dependency between features in order to anticipate counterattacks. One possible avenue that future research could take with the CIC-Darknet2020 dataset is to develop an obfuscation method to exploit RF feature importance, or the weights of a linear SVM. This might better correlate the relationship between classifier response and dataset statistics. We only tested our obfuscation method using our best-performing classifier. It would be interesting to explore how other classifiers respond to similar obfuscation techniques.


Experience report

Fourteen years of R/qtl: Just barely sustainable

Karl W. Broman

Department of Biostatistics & Medical Informatics

University of Wisconsin–Madison

Address for correspondence:

Karl W Broman
---
Department of Biostatistics & Medical Informatics
University of Wisconsin-Madison
2126 Genetics/Biotechnology Center
425 Henry Mall
Madison, Wisconsin 53706
Phone: | 608-262-4633
Email: |<EMAIL_ADDRESS>

## Abstract

R/qtl is an R package for mapping quantitative trait loci (genetic loci that
contribute to variation in quantitative traits) in experimental crosses. Its
development began in 2000. There have been 38 software releases since 2001.
The latest release contains 35k lines of R code and 24k lines of C code, plus
15k lines of code for the documentation. Challenges in the development and
maintenance of the software are discussed. A key to the success of R/qtl is
that it remains a central tool for the chief developer’s own research work,
and so its maintenance is of selfish importance.

## Introduction

If two inbred strains (for example, of plants) show consistent differences in
a quantitative trait, one can be confident that the difference is genetic. An
experimental cross between the two strains can be used to identify the genetic
loci (called quantitative trait loci, QTL) that contribute to the trait
difference: we seek genomic regions for which genotypes are associated with
the trait.

As an illustration, Figure 1 displays the results of QTL analysis with data
from [1], on gravitropism in Arabidopsis. LOD curves measuring the strength of
association between phenotype and genotype are displayed in Figure 1A; a plot
of phenotype vs. genotype at the marker exhibiting the strongest association
is shown in Figure 1B.

Figure 1: Typical analysis results from R/qtl. A: LOD curves across the
genome, measuring association between phenotype and genotype across, and B:
Association between genotype and phenotype at the marker with the strongest
association. The data are from [1]; panel B was created using the R package
beeswarm.

Numerous software packages for QTL analysis are available, some commercial
(e.g., MapQTL and MultiQTL) and some free (e.g. Mapmaker/QTL, and QTL
Cartographer).

My own QTL mapping software, R/qtl [2], is developed as an add-on package to
the widely used general statistical software, R [3]. The software is open
source, licensed under GPL3, and currently hosted on GitHub.

## History

I became interested in QTL mapping in graduate school, twenty years ago. Mark
Neff introduced me to Lander and Botstein’s paper on interval mapping [4],
which remains the most commonly used QTL analysis method.

In the fall of 1999, I joined the Department of Biostatistics at Johns Hopkins
University as an assistant professor. In February, 2000, Gary Churchill
visited me from the Jackson Laboratory, and we discussed our shared interest
in QTL analysis methods and the need for more advanced software. Gary
suggested that we write our own QTL mapping package. He was thinking Matlab,
but I was keen to use R. (R version 1.0 was released the following week.) R
won out over Matlab largely because I developed a working prototype more
quickly. I had recently written some C code implementing the hidden Markov
model (HMM) technology [5] for the treatment of missing genotype information
for QTL analysis in experimental crosses. This served as the starting point
for the package.

Our main goal was for the software to enable the exploration of multiple-QTL
models. We also wanted it to be easily extendible as new methods were
developed. My initial concept was to implement all QTL mapping methods, good
and bad, so that their performance could be compared within a single package.

Much of the development of R/qtl occurred during a three-month sabbatical I
spent at the Jackson Laboratory in Fall, 2001. (It is easy to remember the
year, because I was there on 11 September 2001.) Hao Wu, a software engineer
working with Gary from 2001–2005, contributed a great deal to the core of
R/qtl. In 2009–2010, Danny Arends (with some assistance from Pjotr Prins and
me) incorporated Ritsert Jansen’s MQM code [6, 7, 8], previously available
only in commercial software, as part of R/qtl [9].

## About me

I am an applied statistician. My primary interest is in helping other
scientists answer questions using data. But that generally requires the
development of new statistical methods, and such methods must be implemented
in software. Thus, I spend a considerable amount of time programming.

I have little formal training in programming or software engineering, and I am
not a specialist in computational statistics. But I think it’s important for
an applied statistician to be self-sufficient: We can’t rely on others to
develop the tools we need but must be able to do that ourselves.

## Strengths

R/qtl has a number of strengths. It is comprehensive: It includes
implementations of many QTL mapping methods, and it has a number of tools for
the fit and exploration of multiple-QTL models. It has extensive facilities
for data diagnostics and visualizations. It can be extended, as the results of
important intermediate calculations, that serve as the basis for any QTL
mapping method, are exposed to the user.

The central calculations are coded in C, for speed, but R is used for
manipulating data objects and for graphics. Figure 2 displays the growth of
code in R/qtl over time, as the number lines of code in R and C, as well as in
the R documentation files. Ignoring the documentation files, the code is about
60% R and 40% C.

Figure 2: Numbers of lines of code in released versions of R/qtl over time.
(Yellow and blue correspond to R and C code, respectively; pink is for the R
documentation files.)

Developing the software as a package for R has a number of advantages for the
developer, particularly to make use of R’s facilities for graphics and
input/output, as well as its extensive math/stat software library. The R
packaging system and the Comprehensive R Archive Network (CRAN,
http://cran.r-project.org) simplifies the software installation process and
makes it easy to provide documentation.

The user also benefits by having the QTL mapping software embedded within
general statistical software: for further data analysis before and after the
QTL analysis, and for developing specially-tailored visualizations. R also
provides an excellent interactive environment for data analysis.

A number of related packages are coordinated with R/qtl, using a common data
structure or input file format and some shared functions. Examples include
qtlbim [10], wgaim [11], and dlmap [12].

## Weaknesses

R/qtl also has a number of weaknesses. There has largely been one developer
(who doubles as the support staff), and so many desired features have not been
implemented. We never wrote formal specifications for the internal data
formats, which makes it more difficult for others to contribute code to the
project.

There have been a number of memory management issues, particularly regarding
the copying of large datasets in memory as it is moved from R to C. This can
reduce performance in extremely high-dimensional applications.

The biggest weakness is that the central data structure is too restrictive.
This makes it difficult to extend the software beyond the simplest types of
experimental crosses.

## Some really bad code

More embarrassing than the above weaknesses is that, while R/qtl contains some
good code (like the HMM implementation), there is also some really terrible
code. I’ve learned a great deal about programming while developing R/qtl, but
it’s hard to go back and replace old code. (And some of this code was
obviously bad when written and just should have been constructed with greater
care.)

The worst piece of code (now fixed: see the commit in the git repository on
GitHub, https://github.com/kbroman/qtl/commit/4cd486) involved a really stupid
approach to find the first non-blank element in a vector of character strings.
The code worked fine in small datasets, and so it took a while to discover the
problem. (Open source means everyone can see my stupid mistakes. Version
control means everyone can see every stupid mistake I’ve ever made.)

In many cases, functions aren’t split up into short reusable functions the way
they should be, and there is a lot of repeated code. And some of the functions
have very long lists of arguments, which can be really intimidating to users.

In lots of cases, bugs were fixed by adding little band-aids of code.
Functions get longer and longer to handle more and more special cases, to the
point that they are nearly unreadable. These very long, complex functions are
a barrier to the further development of the software.

The worst offender is scantwo(), for performing a two-dimensional genome scan
for pairs of QTL. The R function is 1402 lines long (that is 4% of the R code
in the package). And the R code is just moving data and results around. The
actual calculations are performed in C, in a series of files that comprise
4725 lines, which is 20% of the C code in the package.

## Version control

In the first eight years of developing R/qtl, I didn’t use formal version
control. Initially, I was simply editing the code in place and saving copies
of releases when I sent them to the R Archive. Incorporating others’
contributions was often a hassle.

In the fall of 2006, Śaunak Sen and I used a Subversion repository to
facilitate the development of our book about QTL mapping and R/qtl [13], but
it wasn’t until January, 2008, that I began to use subversion for R/qtl, as
well. In February, 2009, Pjotr Prins helped me to switch to git and to place
the R/qtl repository on GitHub.

I’m not sure how I managed the project before adopting version control.
Version control makes it so easy to try things out without fear of breaking
working code. And collaboration on software development is terribly tedious
and inefficient without version control.

## User support

Helping users to solve their research problems can be quite rewarding. But it
can also require considerable patience. I respond to several queries per week
about R/qtl, mainly through a Google Group that was started in September,
2004, but many others directly via email.

I had hoped that the R/qtl discussion group would become a network of folks
answering each others’ questions about R/qtl and QTL mapping, but I am
basically the only one answering questions. I think that’s partly because, to
avoid spam, I moderate all messages, at which point I generally try to answer
them. Thus, participants see my answer at basically the same time as they see
the original question. But at least these electronic conversations are public
and searchable.

I’ve helped many people revise their data files for use with R/qtl. It’s very
hard to predict what might go wrong in data import and to provide appropriate
error messages.

Questions are often frustrating. They may provide too little detail (“I tried
such-and-such but _it didn’t work_.”) or too much (“Could you please look at
the attached 25-page Word document containing code and output and tell me if
I’m doing something wrong?”) Many questions are not so much about the software
but are more general scientific or data analysis questions.

Many questions would be answered by a careful reading of the documentation,
but software documentation is often dreadfully boring. I’ve learned that the
most popular documentation are the short tutorials with practical examples
clearly illustrating important tools. These are time-consuming to create (and
maintain).

Some things are _possible_ in R/qtl but are not well documented and not really
feasible except for users with considerable programming and analysis
experience. It feels wrong to say, “That _is_ possible, but I don’t have the
time to explain how.”

It was helpful to have written a book about QTL mapping and R/qtl [13], but
it’s frustrating to watch the publisher nearly double its street price. If I
could do it again, I would self-publish.

## Sustainable academic software

Why put so much effort into software development? The principal advantage is
the software itself: QTL analysis is the focus of much of my own research, and
software that is easy for others to use is also easy for me to use. Second,
R/qtl provides a platform for me to distribute implementations of QTL mapping
methods that I develop, such as the two-part model of [14] and the model
selection methods of [15] and [16]. Third, the software has led to many
interesting consulting-type interactions and a good number of collaborations
(leading to co-authorship on papers and some grant support). Fourth, my own
methods grant is much more attractive with a successful software aim. Finally,
the software supports others’ research, and my primary duty as a statistician
and an academic is to help others.

I think the key to the sustainability of piece of scientific software is that
the developer continues to use the software for his/her own research. So
often, the developer moves on to other research problems, and the software is
orphaned.

But software development requires considerable time and support. Often,
researchers don’t recognize the many indirect benefits of building useful
software tools, nor the long-term benefits of devoting extra time to the
construction of _good_ software. The reviewers of grant proposals may
recognize the value of software development, but generally it can not be the
sole research effort, and funding agencies often over-emphasize innovation,
which can make it difficult to obtain support for software maintenance.
Moreover, academic culture encourages researchers to create independent
software packages rather than contribute to existing ones.

The biggest open question, for me personally, is how to encourage the
formation of a community around a given software development effort: for
coding, for answering users’ questions, and for developing improved
documentation and tutorials.

## Future

As Fred Brooks said in _The Mythical Man Month_ , “Plan to throw one away; you
will anyhow.” [17], Ch. 11) This was restated by Eric Raymond in _The
Cathedral and The Bazaar_ : “You often don’t really understand the problem
until after the first time you implement a solution. The second time, maybe
you know enough to do it right.” [18]

In collaboration with Danny Arends and Pjotr Prins, I’ve initiated a
reimplementation of R/qtl, with a focus on high-dimensional data and more
modern cross designs, such as the Collaborative Cross [19] and the Mouse
Diversity Outbred Population [20]. We will also focus on interactive graphical
tools, implemented with the Javascript library, D3.

There has been a renewal of interest in QTL analysis, particularly with the
growth of eQTL analysis, in which genome-wide gene expression measures are
treated as phenotypes (see, for example, [21, 22]). We hope that our new
software will become a popular platform for the analysis of large-scale eQTL
data, as R/qtl has been for more traditional-sized QTL projects.

## Acknowledgments

This work was supported in part by NIH grant R01 GM074244. Christina
Kendziorski generously provided comments to improve the manuscript. Numerous
people have contributed to R/qtl over the years: Danny Arends, Gary Churchill,
Ritsert Jansen, Steffen Moeller, Pjotr Prins, Śaunak Sen, Laura Shannon, Hao
Wu, and Brian Yandell.

## References

  * [1] Moore CR, Johnson LS, Kwak IY, Livny M, Broman KW, Spalding EP.  High-Throughput Computer Vision Introduces the Time Axis to a Quantitative Trait Map of a Plant Growth Response.  Genetics. 2013;195:1077–1086.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3813838/.
  * [2] Broman KW, Wu H, Sen S, Churchill GA.  R/qtl: QTL mapping in experimental crosses.  Bioinformatics. 2003;19:889–890.  Available from: http://bioinformatics.oxfordjournals.org/content/19/7/889.long.
  * [3] R Core Team.  R: A language and environment for statistical computing; 2013.  Available from: http://www.r-project.org.
  * [4] Lander ES, Botstein D.  Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps.  Genetics. 1989;121:185–199.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1203601/.
  * [5] Baum LE, Petrie T, Soules G, Weiss N.  A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains.  Ann Math Statist. 1970;41:164–171.  Available from: http://www.jstor.org/stable/2239727.
  * [6] Jansen RC.  Interval mapping of multiple quantitative trait loci.  Genetics. 1993;135:205–211.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1205619/.
  * [7] Jansen RC.  Controlling the type I and type II errors in mapping quantitative trait loci.  Genetics. 1994;138:871–881.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1206235/.
  * [8] Jansen RC, Stam P.  High resolution of quantitative traits into multiple loci via interval mapping.  Genetics. 1994;136:1447–1455.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1205923/.
  * [9] Arends D, Prins P, Jansen RC, Broman KW.  R/qtl: high throughput multiple QTL mapping.  Bioinformatics. 2010;26:2990–2992.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2982156/.
  * [10] Yandell BS, Mehta T, Banerjee S, Shriner D, Venkataraman R, Moon JY, et al.  R/qtlbim: QTL with Bayesian interval mapping in experimental crosses.  Bioinformatics. 2007;23:641–643.  Available from: http://bioinformatics.oxfordjournals.org/content/23/5/641.long.
  * [11] Taylor J, Verbyla A.  R package wgaim: QTL analysis in bi-parental populations using linear mixed models.  J Statist Soft. 2011;40:1–18.  Available from: http://www.jstatsoft.org/v40/i07/paper.
  * [12] Huang BE, Shah R, George AW.  dlmap: An R package for mixed model QTL and association analysis.  J Statist Soft. 2012;50:1–22.  Available from: http://www.jstatsoft.org/v50/i06/paper.
  * [13] Broman KW, Sen S.  A guide to QTL mapping with R/qtl.  Springer; 2009.
  * [14] Broman KW.  Mapping quantitative trait loci in the case of a spike in the phenotype distribution.  Genetics. 2003;163:1169–1175.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1462498/.
  * [15] Broman KW, Speed TP.  A model selection approach for the identification of quantitative trait loci in experimental crosses.  J R Statist Soc B. 2002;64:641–656.  Available from: http://www.jstor.org/stable/3088807.
  * [16] Manichaikul A, Moon JY, Sen S, Yandell BS, Broman KW.  A model selection approach for the identification of quantitative trait loci in experimental crosses, allowing epistasis.  Genetics. 2009;181:1077–86.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2651044/.
  * [17] Brooks FP.  Mythical Man Month: Essays on Software Engineering.  Addison-Wesley; 1975.
  * [18] Raymond ES.  The Cathedral & the Bazaar.  O’Reilly, Sebastopol, CA; 1999.  Available from: http://www.catb.org/~esr/writings/cathedral-bazaar/cathedral-bazaar/.
  * [19] The Complex Trait Consortium.  The Collaborative Cross, a community resource for the genetic analysis of complex traits.  Nat Genet. 2004;36:1133–1137.  Available from: http://www.nature.com/ng/journal/v36/n11/full/ng1104-1133.html.
  * [20] Svenson KL, Gatti DM, Valdar W, Welsh CE, Cheng R, Chesler EJ, et al.  High-resolution genetic mapping using the Mouse Diversity outbred population.  Genetics. 2012;190:437–447.  Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3276626/.
  * [21] Jansen RC, Nap JP.  Genetical genomics: the added value from segregation.  Trends Genet. 2001;17:388–391.  Available from: http://www.cell.com/trends/genetics/abstract/S0168-9525(01)02310-1.
  * [22] Tesson BM, Jansen RC.  eQTL analysis in mice and rats.  In: DiPetrillo K, editor. Cardiovascular Genomics, Methods in Molecular Biology 573. 10. Humana Press; 2009. p. 285–309.