Dataset Exploration and Analysis

The cover image is generated by ChatGPT, depicting data analysis. The prompt used is “A modern data analysis concept illustration featuring a futuristic workspace. The image includes multiple data charts, graphs, and dashboards displayed on transparent holographic screens. A diverse team of analysts and data scientists collaborate, analyzing trends and insights on large monitors. The scene has a sleek, high-tech atmosphere with glowing blue and purple hues, reflecting a professional and cutting-edge environment.” 。

Introduction

The dataset used in this analysis is the Swiss bank notes dataset, which is utilized by banks to distinguish between genuine and counterfeit old Swiss Francs. This study downloads the bank2.dat dataset from https://github.com/QuantLet/MVA/tree/master/QID-1530-MVAscabank56 and conducts the analysis using the R programming language.

Dataset

From the textbook, we know that the first 100 records in this dataset correspond to genuine banknotes, while the last 100 records correspond to counterfeit banknotes. The dataset consists of 200 observations and 6 variables.

Variable Description

According to Table 22.3, the dataset contains 6 variables, described as follows:

We read the dataset in R and create a new variable, genuine, to indicate whether the banknote is genuine (1 for genuine, 0 for counterfeit). The corresponding R code is as follows:

1
2
3
4
5

data = read.table('bank2.dat')
colnames(data) = c('X1','X2','X3','X4','X5','X6')
data = data.frame(data, genuine = c(rep(1, 100), rep(0, 100)))
attach(data)
head(data)

The Dataset

A portion of the dataset is shown below.

Summary

We can use the following code to quickly inspect the dataset.

1
2
3
4
5

print('Quartiles:')
summary(data)

print('NAs:')
sapply(data, function(x) sum(is.na(x)))

A brief summary of the dataset is as follows.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

[1] "Quartiles:"
       X1              X2              X3              X4               X5              X6           genuine   
 Min.   :213.8   Min.   :129.0   Min.   :129.0   Min.   : 7.200   Min.   : 7.70   Min.   :137.8   Min.   :0.0  
 1st Qu.:214.6   1st Qu.:129.9   1st Qu.:129.7   1st Qu.: 8.200   1st Qu.:10.10   1st Qu.:139.5   1st Qu.:0.0  
 Median :214.9   Median :130.2   Median :130.0   Median : 9.100   Median :10.60   Median :140.4   Median :0.5  
 Mean   :214.9   Mean   :130.1   Mean   :130.0   Mean   : 9.418   Mean   :10.65   Mean   :140.5   Mean   :0.5  
 3rd Qu.:215.1   3rd Qu.:130.4   3rd Qu.:130.2   3rd Qu.:10.600   3rd Qu.:11.20   3rd Qu.:141.5   3rd Qu.:1.0  
 Max.   :216.3   Max.   :131.0   Max.   :131.1   Max.   :12.700   Max.   :12.30   Max.   :142.4   Max.   :1.0  

 [1] "NAs:"
     X1      X2      X3      X4      X5      X6 genuine 
      0       0       0       0       0       0       0 

From the above response, we can see that there are no missing values in this dataset.

Correlation Coefficient

The correlation coefficient is a method used to measure the linear relationship between two variables. It is calculated as follows. $$ \rho_{X, Y} = \frac{cov(X, Y)}{\sigma_X \sigma_Y} = \frac{E \left[ (X - \mu_X )(Y - \mu_Y) \right]}{\sigma_X \sigma_Y}, $$ Where $cov(X, Y)$ is the covariance between $X$ and $Y$, $\mu_X$ and $\mu_Y$ are the population means of $X$ and $Y$, and $\sigma_X$ and $\sigma_Y$ are the population standard deviations of $X$ and $Y$.

The above is known as the population correlation coefficient. In a dataset, by estimating the sample covariance and standard deviation, we can calculate the sample correlation coefficient, which is given by the following formula. $$ r = \frac{\sum_{i=1}^n(X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum_{i=1}^n(X_i - \bar{X})^2} \sqrt{\sum_{i=1}^n(Y_i - \bar{Y})^2}} = \frac{1}{n - 1}\sum_{i=1}^n \left(\frac{(X_i - \bar{X})}{\sigma_X}\right) \left(\frac{(Y_i - \bar{Y})}{\sigma_Y}\right), $$ Where $X_i$ and $Y_i$ represent the values of the $i$-th sample in the dataset, and $\bar{X}$ and $\bar{Y}$ are the sample means of $X$ and $Y$, respectively.

In R, the cor(dataset) function can be used to generate the correlation coefficients between each variable. Additionally, various packages can be used to visualize the correlation coefficients. Below, we use the ggcorrplot package to visualize the correlation matrix. The correlation coefficients are displayed, with blue representing negative correlations and red representing positive correlations. The size of the points increases with the strength of the correlation, while points with correlations close to 0 become smaller.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

library(ggcorrplot)

cor_matrix = cor(data)

ggcorrplot(cor_matrix, 
           method = "circle",
           type = "full",
           lab = TRUE,
           lab_size = 3,
           colors = c("blue", "white", "red"),
           outline.color = "gray",
           legend.title = "Correlation",
           show.legend = TRUE)

From the above plot, we can observe that $X_6$ has a strong positive correlation in distinguishing between genuine and counterfeit banknotes, while $X_1$ shows little to no linear correlation.

Scatter Plot

A scatter plot is one of the simplest ways to visualize data. You can easily create a scatter plot using plot(x_coordinates, y_coordinates). If there are many variables, you can use pairs to create scatter plots for every pair of variables. Here’s how to do it.

1

pairs(data, upper.panel = NULL)

From the above plot, we can observe that there seems to be a linear relationship between $X_2$, $X_3$, $X_4$, $X_5$, and $X_6$, while $X_1$ does not show a strong relationship with the other variables. Also, since genuine is a categorical variable, its interactions with other variables are always confined to 0 and 1.

2D Scatter Plot

Of course, we can also plot a scatter plot for each pair of variables individually. Below, we will plot the scatter plot for each pair of data and use the genuine variable to distinguish between genuine and counterfeit notes with different colors.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

pairs_plot = function(x, y, x_name, y_name){
  par(mar = c(5, 5, 4, 10))
  plot(x, y, main = "Swiss bank notes", xlab = x_name, ylab = y_name)
  points(x[1:100], y[1:100], col = "blue")
  points(x[101:200], y[101:200], col = "red")
  legend("topright", legend = c("Genuine", "Counterfeit"), col = c("blue", "red"), pch = 1, xpd = TRUE, inset = c(-0.3, 0))
}

num = combn(1:6, 2)
for (i in 1:ncol(num)) {
  x = paste0("X", num[1, ][i])
  y = paste0("X", num[2, ][i])
  pairs_plot(get(x), get(y), x, y)
}

We can obtain the following scatter plots.

From the above pairwise scatter plots, we can observe that some plots, such as those of $X_1$ and $X_2$, $X_1$ and $X_3$, $X_1$ and $X_5$, and $X_2$ and $X_5$, show very fuzzy boundaries and do not provide a clear distinction between genuine and counterfeit notes. On the other hand, plots such as $X_2$ and $X_6$, $X_3$ and $X_6$, $X_4$ and $X_5$, $X_4$ and $X_6$, and $X_5$ and $X_6$ demonstrate more distinct separations. Using these pairwise variables, it is easier to distinguish between genuine and counterfeit notes.

With these more distinguishable scatter plots, methods such as Decision Trees, Support Vector Machines (SVM), and K-Nearest Neighbors (KNN) can be used to classify and predict whether new banknotes (data) are genuine or counterfeit.

3D Scatter Plot

We can also create a 3D scatter plot. Here is the code to do so.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

library(lattice)
library(grid)

groups = data[, 7]

pch_types = c(1, 1)              
colors = c("blue", "red")        
pch_group = pch_types[ifelse(groups == 0, 1, 2)]
col_group = colors[ifelse(groups == 0, 1, 2)]

combos = combn(1:6, 3)

n_rows = 1        
n_cols = 1        
plots_per_page = n_rows * n_cols   

for (page in seq(1, ncol(combos), by = plots_per_page)) {
  grid.newpage()

  for (i in 0:(plots_per_page - 1)) {
    index = page + i
    if (index > ncol(combos)) break

    x_name = paste0("X", combos[1, index])
    y_name = paste0("X", combos[2, index])
    z_name = paste0("X", combos[3, index])

    x_data = data[[x_name]]
    y_data = data[[y_name]]
    z_data = data[[z_name]]

    p = cloud(z_data ~ x_data * y_data,
          pch = pch_group,
          col = col_group,
          cex = 1.2,
          ticktype = "detailed",
          main = paste("Swiss bank notes:", x_name, y_name, z_name),
          screen = list(z = -90, x = -90, y = 45),
          scales = list(arrows = FALSE, col = "black", distance = 1, cex = 0.5),
          xlab = list(x_name, rot = -10, cex = 1.2),
          ylab = list(y_name, rot = 10, cex = 1.2),
          zlab = list(z_name, rot = 90, cex = 1.1),

          par.settings = list(
            axis.line = list(col = "black"),    
            box.3d = list(
              col = c(
                "black", "transparent", "transparent","black", "black", "transparent", "transparent", "transparent",   
                "transparent", "transparent", "transparent", "transparent"))),

          key = list(
            space = "right",
            points = list(pch = pch_types, col = colors, cex = 1.5),
            text = list(c("Genuine", "Counterfeit")),
            border = FALSE
          )
        )

    print(p, split = c((i %% n_cols) + 1, (i %/% n_cols) + 1, n_cols, n_rows), more = TRUE)
  }
}

From the 3D scatter plot, we can better observe the interactions between any three sets of data.

We can also use the plotly package to create an interactive 3D scatter plot.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

library(plotly)
library(htmlwidgets)

x_data = data$X1
y_data = data$X2
z_data = data$X6
groups = data[, 7]

groups = factor(groups, levels = c(0, 1), labels = c("Genuine", "Counterfeit"))

x_range <- range(x_data)
y_range <- range(y_data)
z_range <- range(z_data)

x_scaled = (x_data - x_range[1]) / (x_range[2] - x_range[1])
y_scaled = (y_data - y_range[1]) / (y_range[2] - y_range[1])
z_scaled = (z_data - z_range[1]) / (z_range[2] - z_range[1])

tick_positions = seq(0, 1, length.out = 5)
tick_labels_x = round(seq(x_range[1], x_range[2], length.out = 5), 1)
tick_labels_y = round(seq(y_range[1], y_range[2], length.out = 5), 1)
tick_labels_z = round(seq(z_range[1], z_range[2], length.out = 5), 1)

hover_text <- paste0(
  "X1: ", round(x_data, 1), "<br>",
  "X2: ", round(y_data, 1), "<br>",
  "X6: ", round(z_data, 1)
)

fig <- plot_ly(
  x = x_scaled,
  y = y_scaled,
  z = z_scaled,
  color = groups,
  colors = c("blue", "red"),
  type = "scatter3d",
  mode = "markers",
  marker = list(size = 5),
  text = hover_text,
  hoverinfo = "text" 
)

fig <- fig %>% layout(
  title = "3D Scatter Plot of X1, X2, X6",
  scene = list(
    xaxis = list(title = "X1", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_x),
    yaxis = list(title = "X2", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_y),
    zaxis = list(title = "X6", range = c(0, 1), tickvals = tick_positions, ticktext = tick_labels_z)
  )
)

saveWidget(fig, "3d_scatter_plot.html")

If you are unable to view the interactive 3D scatter plot or need to view it in full screen, please click here to access it.

Box Plot

A box plot, also known as a box-and-whisker plot, arranges all values in the dataset from smallest to largest and divides them into four equal parts. The values at the three dividing points are the quartiles. From a box plot, we can easily identify the extreme values and the quartiles. We can use the following code to plot a box plot.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

name = NULL
for (i in 1:6) {
  name = c(name, paste0("X", i, "_Genuine"))
  name = c(name, paste0("X", i, "_Counterfeit"))
}

groups = list()
for (i in 1:length(name)){
  if (i %% 2 == 1) {
    groups[[name[i]]] = data[1:100, as.integer(i / 2) + 1]
  } else {
    groups[[name[i]]] = data[101:200, as.integer(i / 2)]
  }
}

means <- sapply(groups, mean)

par(mfrow=c(1,3))
for (i in seq(1, length(name) ,by = 2)){
  boxplot(groups[i:(i + 1)], names = name[i:(i + 1)], frame = TRUE, main = "Swiss bank notes", cex.axis=0.99)
  
  for (j in 0:1) {
    lines(c(j + 0.6, j + 1.4), rep(means[[name[(i + j)]]], 2), lty = "dotted", lwd = 1.2)
  }
}

In each of the plots, the dashed line represents the mean of the data, while the solid line represents the second quartile (median, Q2). The upper and lower limits of the box represent the third quartile (Q3) and the first quartile (Q1), respectively. The upper and lower boundaries show the maximum and minimum values, excluding extreme values (outliers). The interquartile range (IQR) is the difference between Q3 and Q1. If a data point is greater or less than 1.5 times the IQR from the box, it is considered an outlier.

Upon observing the box plots, we can see that for $X_1$, $X_2$, $X_3$, and $X_5$, the box plots for genuine and counterfeit notes do not show significant differences, making these variables less suitable for standalone analysis. Although the box plot for $X_4$ indicates a larger difference in mean and median between the two categories, there is still a substantial overlap in the box plot for genuine notes, which makes it less suitable for standalone analysis as well. However, for $X_6$, the box plot clearly shows a distinct difference between genuine and counterfeit notes, with almost no overlap in the data range, making it suitable for analysis and differentiation between genuine and counterfeit notes.

Histogram

A histogram is a two-dimensional statistical chart that groups data into intervals and represents the number of data points within each interval with the height of bars. It provides an intuitive way to observe the distribution of data. In R, we can use the hist() function to plot histograms. Below, we plot histograms for each variable and distinguish between genuine and counterfeit notes using blue and red colors.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

data_Genuine = data[1:100, 1:6]   
data_Counterfeit = data[101:200, 1:6]

for (i in 1:6) {
  breaks_seq = seq(min(data[, i]), max(data[, i]), by = 0.1)
  min_x = floor(min(data[, i]))
  max_x = ceiling(max(data[, i]))

  hist(data_Genuine[, i],
       breaks = breaks_seq,
       col = rgb(0, 0, 1, 0.5),   
       border = "blue",
       xlim = c(min_x, max_x),
       ylim = c(0, max(table(cut(data[, i], breaks_seq))) + 5),
       main = "Histogram",
       xlab = paste0("X", i),
       axes = FALSE) 
  
  hist(data_Counterfeit[, i],
       breaks = breaks_seq,
       col = rgb(1, 0, 0, 0.5),   
       border = "red",
       add = TRUE)
  
  axis(side = 1, at = seq(min_x, max_x, by = 1))
  axis(side = 2)
  
  legend("topright",
         legend = c("Genuine", "Counterfeit", "Overlap"),
         fill = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5), rgb(0.6, 0, 0.4, 0.7)),
         border = c("blue", "red", "purple"))
}

From the histograms, we can observe that in the variable $X_6$, the overlapping part is the smallest, accounting for only about 2% of the total. In contrast, the overlapping part in $X_1$ is the largest, accounting for approximately 73% of the total. This indicates that, compared to $X_1$, using $X_6$ or $X_4$ for classification could provide higher discriminative power. The overlap in $X_6$ is only 2%, and in $X_4$, it is only 15%, meaning these variables can more effectively distinguish between genuine and counterfeit notes. On the other hand, the overlap in $X_1$ reaches 76%, suggesting that $X_1$ performs weaker in differentiating between the two categories. Therefore, choosing $X_6$ may lead to better recognition of genuine and counterfeit notes.

Kernel Density Estimation (KDE)

Kernel density estimation is a method used to estimate the probability density function (PDF) of a random variable. For each data point in a dataset, KDE draws a local distribution based on a kernel function, and then sums these local distributions to form the overall estimated density function. Compared to histograms, kernel density estimation provides a smoother and continuous distribution curve.

Assuming $x_1, x_2, \cdots, x_n$ are independent and identically distributed (i.i.d.) samples, we can estimate the kernel density using the following formula:

$$ \hat{f}_h (x) = \frac{1}{n} \sum_{i=1}^{n} K_h (x - x_i) = \frac{1}{nh} \sum_{i=1}^{n} K \left( \frac{x - x_i}{h} \right), $$

Here, $K$ is a non-negative kernel function, and $h > 0$ is the smoothing parameter (bandwidth).

In this case, we use the Biweight kernel function, which is given by the following formula: $$ K(u) = \begin{cases} \frac{16}{15} (1 - u^2)^2, & \text{for } |u| \leq 1; \\ 0, & \text{for } |u| > 1. \end{cases} $$

We can use the following code to plot the estimated probability density function.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

library(KernSmooth)

for (i in 1:6) {
  
  data_Genuine = data[1:100, i]   
  data_Counterfeit = data[101:200, i]  

  fh1 = bkde(data_Genuine, kernel = "biweight")  
  fh2 = bkde(data_Counterfeit, kernel = "biweight")  

  x_min = floor(min(c(fh1$x, fh2$x)))
  x_max = ceiling(max(c(fh1$x, fh2$x)))
  y_max = max(c(fh1$y, fh2$y)) * 1.1

  x_ticks <- seq(x_min, x_max, by = 1)

  plot(fh1, type = "l", lwd = 2,
       xlab = "Counterfeit / Genuine", 
       ylab = paste0("Density estimates for X", i, sep = ""),
       col = "blue", main = paste("Swiss bank notes - X", i),
       xlim = c(x_min, x_max), ylim = c(0, y_max),
       axes = FALSE) 
  
  lines(fh2, lty = "dotted", lwd = 2, col = "red")

  axis(side = 1, at = x_ticks, labels = x_ticks)
  axis(side = 2)
  box()

  legend("topright",
         legend = c("Genuine", "Counterfeit"),
         col = c("blue", "red"),
         lty = c("solid", "dotted"),
         lwd = 2,
         cex = 0.8)
}

Compared to histograms, kernel density estimation can approximate the probability density function in a smoother, continuous way. It is important to note that for both kernel density estimation and histograms, choosing an appropriate bandwidth is a crucial factor that affects the graphical result. Too much or too little bandwidth can lead to overly fragmented or coarse graphs, making it harder to discern the overall trend.

The kernel density estimation plots created using $X_1$ to $X_6$ show that the degree of overlap between the blue and red curves reflects the variable’s ability to discriminate between classes. The more clearly the curves separate, the better the variable is at distinguishing between true and counterfeit notes. From the plots, we can draw the following conclusions:

• $X_1$: The peak is slightly skewed, but the two curves almost overlap, with a large overlapping range, indicating weak discriminative ability.
• $X_2$: The peak is slightly skewed, with a wide overlapping range. The main peaks differ but the difference is limited, suggesting weak discriminative ability.
• $X_3$: There is a noticeable difference in the distribution curves for true and counterfeit notes, with distinct peaks, indicating moderate discriminative ability.
• $X_4$: The distribution ranges for true and counterfeit notes are separate with minimal overlap, indicating excellent discriminative ability.
• $X_5$: The distribution shapes and ranges for true and counterfeit notes are similar, but with different peaks, indicating moderate discriminative ability.
• $X_6$: The distribution ranges for true and counterfeit notes are separate with very little overlap, indicating excellent discriminative ability.

Andrews’ Curves

Andrews’ Curves utilize the expansion form of the Fourier series for visualizing multidimensional data, transforming each data point into a continuous curve. This allows for an intuitive representation of similarities and differences between samples in a two-dimensional coordinate system.

Each multivariate observation $X_i = \begin{pmatrix} X_{i,1} & X_{i,2} & \cdots & X_{i,p} \end{pmatrix}$ can be converted into a parameterized curve $f_i(t)$. Depending on whether the dimensionality $p$ is odd or even, its expansion form is given by: $$ f_i(t) = \begin{cases} \frac{X_{i,1}}{\sqrt{2}} + X_{i,2} \sin(t) + X_{i,3} \cos(t) + \cdots + X_{i,p-1} \sin \left( \frac{p-1}{2} t \right) + X_{i,p} \cos \left( \frac{p-1}{2} t \right), & \text{for } p \text{ odd}; \\ \frac{X_{i,1}}{\sqrt{2}} + X_{i,2} \sin(t) + X_{i,3} \cos(t) + \cdots + X_{i,p} \sin \left( \frac{p}{2} t \right), & \text{for } p \text{ even}. \end{cases} $$

We can construct the following code to plot Andrews’ Curves.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

library(tourr)

x = data[1:200, ]
y = NULL

for (i in 1:6) {
  z = (x[, i] - min(x[, i])) / (max(x[, i]) - min(x[, i]))  # zero-one scaling
  y = cbind(y, z)
}

Type = data[, 7]
f = as.integer(Type)
grid = seq(0, 2 * pi, length = 1000)

plot(grid, 
     andrews(y[1, ])(grid), 
     type = "l", 
     lwd = 1.2, 
     main = "Andrews' curves (All Bank data)", 
     axes = FALSE, 
     frame = TRUE, 
     ylim = c(-0.5, 0.6), 
     ylab = "", 
     xlab = "", 
     col = ifelse(f[1] == 1, "blue", "red"), 
     lty = ifelse(f[1] == 1, "solid", "dotted"))

for (i in 2:200) {
  lines(grid, 
        andrews(y[i, ])(grid), 
        col = ifelse(f[i] == 1, "blue", "red"), 
        lwd = 1.2, 
        lty = ifelse(f[i] == 1, "solid", "dotted"))
}

axis(side = 2, at = seq(-0.5, 0.6, 0.2), labels = seq(-0.5, 0.6, 0.2))
axis(side = 1, at = seq(0, 7, 1), labels = seq(0, 7, 1))

legend("topright", 
       legend = c("Genuine", "Counterfeit"), 
       col = c("blue", "red"), 
       lwd = 1.5, 
       lty = c("solid", "dotted"))

From the Andrews’ Curves shown in the figure, we can observe that the blue curve (representing genuine banknotes) primarily concentrates in the middle fluctuating region, remaining relatively stable overall. On the other hand, the red curve (representing counterfeit banknotes) exhibits higher variability in different sections, with a wider range of fluctuations. At $t \approx 2$, $t \approx 4$, and $t \approx 6$, the counterfeit curve displays more dispersion and amplitude variation in these regions, indicating that there is greater variability among the samples in this group.

Logistic Regression

For general classification problems, we can use logistic regression and random forests to build models.

Logistic regression is a binary classification model represented by the conditional probability $P(Y|X)$, which follows a parametric logistic distribution. The random variable $X$ is a real number, while $Y$ takes values of 0 or 1.

For the banknote authentication dataset, we use a logistic regression model, setting genuine as the dependent variable and $X_1, \cdots, X_6$ as the independent variables. In this implementation, 70% of the data is used for training, and 30% is used for testing. The analysis is as follows.

Data Analysis and Model Summary

A brief data analysis and model summary are as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# summary
set.seed(123456)
split = sample.split(data,SplitRatio = 0.7)
train_reg = subset(data, split =='TRUE')
test_reg = subset(data, split == 'FALSE')

logistic_model = glm(factor(genuine) ~ X1+X2+X3+X4+X5+X6, data = train_reg,
                     family = binomial(link = 'logit'), 
                     control = list(maxit=1000))

logistic_model
summary(logistic_model)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

Warning:  glm.fit: Fitted probabilities numerically equal to 0 or 1
Call:  glm(formula = factor(genuine) ~ X1 + X2 + X3 + X4 + X5 + X6, 
    family = binomial(link = "logit"), data = train_reg, control = list(maxit = 1000))

Coefficients:
(Intercept)           X1           X2           X3           X4           X5           X6  
  -5619.203       15.488        4.479      -17.391      -15.337      -22.686       30.986  

Degrees of Freedom: 113 Total (i.e. Null);  107 Residual
Null Deviance:	    158 
Residual Deviance: 4.51e-10 	AIC: 14

Call:
glm(formula = factor(genuine) ~ X1 + X2 + X3 + X4 + X5 + X6, 
    family = binomial(link = "logit"), data = train_reg, control = list(maxit = 1000))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.619e+03  2.773e+08       0        1
X1           1.549e+01  2.267e+05       0        1
X2           4.479e+00  8.012e+05       0        1
X3          -1.739e+01  1.287e+06       0        1
X4          -1.534e+01  3.596e+05       0        1
X5          -2.269e+01  1.828e+05       0        1
X6           3.099e+01  1.682e+05       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.5800e+02  on 113  degrees of freedom
Residual deviance: 4.5102e-10  on 107  degrees of freedom
AIC: 14

Number of Fisher Scoring iterations: 27

From the output, we can see that the model’s AIC is 14, and the coefficient descriptions are as follows:

• The coefficients for $X_1$, $X_2$, and $X_6$ are positive, indicating that higher values increase the probability of the banknote being genuine.
• The coefficients for $X_3$, $X_4$, and $X_5$ are negative, indicating that higher values decrease the probability of the banknote being genuine.

Confusion Matrix and ROC Curve

We use the following code to make predictions.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

predict_reg = predict(logistic_model,test_reg,type='response')
predict_reg = ifelse(predict_reg >0.5,1,0)
table(test_reg$genuine,predict_reg)
missing_classerr = mean(predict_reg!=test_reg$genuine)
print(paste('Accuracy=',1-missing_classerr))

ROCPred = prediction(predict_reg, test_reg$genuine)
ROCPer = performance(ROCPred, measure = 'tpr',x.measure = 'fpr')

auc = performance(ROCPred,measure='auc')
auc = auc@y.values[[1]]
auc

plot(ROCPer, colourise=TRUE, print.cuttoffs.at=seq(0.1,by=0.1), col=2, lwd=2, main='ROC Curve')
abline(a = 0, b = 1)
auc = round(auc, 4)
legend(.6, .4, auc, title = 'AUC', cex = 1)

1
2
3
4
5
6

   predict_reg
     0  1
  0 42  0
  1  1 43
[1] "Accuracy= 0.988372093023256"
[1] 0.9886364

The confusion matrix shows that the model has a high classification accuracy for both classes. It achieves a perfect classification rate for counterfeit banknotes (Class 0) and a very low error rate for genuine banknotes (Class 1).

The ROC curve (Receiver Operating Characteristic curve) illustrates the trade-off between benefit (True Positive Rate, TPR) and cost (False Positive Rate, FPR) for the classifier. The AUC (Area Under Curve) represents the area under the ROC curve, serving as a common statistical measure of a classifier’s predictive performance.

The AUC value of 0.9886 confirms that the logistic regression model performs exceptionally well in classifying banknotes on the test set (86 observations). This aligns with the confusion matrix results, which show 98.8% accuracy, TPR = 0.977, and FPR = 0.

Random Forest

Random Forest is a learner composed of many different decision trees. The principle behind it is to combine multiple “weak learners” to construct a stronger model. The basic idea of a “strong learner” is to combine several CART trees and introduce randomly assigned training data to improve the computational results.

Data Analysis and Model Summary

A brief data analysis and model summary are as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

library(randomForest)

# summary
set.seed(12345)
data$genuine = as.factor(data$genuine)
split = sample.split(data,SplitRatio = 0.6)
train_rf = subset(data, split =='TRUE')
test_rf = subset(data, split == 'FALSE')
rf = randomForest(genuine~ X1+X2+X3+X4+X5+X6,data =train_rf,ntree=500)
print(rf)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

Call:
 randomForest(formula = genuine ~ X1 + X2 + X3 + X4 + X5 + X6,      data = train_rf, ntree = 500) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 0.88%
Confusion matrix:
   0  1 class.error
0 56  0  0.00000000
1  1 56  0.01754386

From the model’s performance, we can see that the 0.88% OOB (Out-of-Bag) error rate indicates excellent performance. On average, only 1 out of every 113 predictions is incorrect.

The confusion matrix shows that the model has highly accurate classification for both classes. It achieves a perfect classification rate for Class 0 (0% error rate) and a very low error rate for Class 1 (1.75%).

Important Variables

1
2

importance(rf)
varImpPlot(rf)

1
2
3
4
5
6
7

   MeanDecreaseGini
X1         1.078579
X2         2.833452
X3         4.902606
X4        14.906531
X5         5.751345
X6        26.515655

Based on feature selection, we can see that $X_6$ and $X_4$ have the greatest influence. Other features fall within the range of 0 to 10.

ROC curve

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

pred1 = predict(rf,test_rf, type='prob')
perf = prediction(pred1[,2], test_rf$genuine)
auc = performance(perf, 'auc')
auc_value = auc@y.values[[1]]

pred3 = performance(perf, 'tpr', 'fpr')

plot(pred3, main="ROC Curve for Random Forest", col=2, lwd=2)
abline(a=0, b=1, lwd=2, lty=2, col='gray')

text(0.6, 0.2, paste("AUC =", round(auc_value, 3)), col=2, cex=1.2)

The ROC curve and AUC = 1 indicate that the Random Forest model is a perfect classifier for the dataset, achieving 100% true positive rate and 0% false positive rate.

Environment

• Operating System: Windows 11 24H2
• Programming Language: R 4.4.2

Further Learning

• The R Notebook HTML file used in this article.

References

• Wolfgang Karl Härdle, Léopold Simar. (2015). Applied Multivariate Statistical Analysis (4th ed.). Springer. https://link.springer.com/book/10.1007/978-3-662-45171-7

• Quantnet. (2025, January 27). MVA. Github. https://github.com/QuantLet/MVA