Overview

This Rmarkdown documents contains the code and 2D simulations of the error propagation and classification uncertainty estimation using Monte Carlo simulations for two generic classes with linear separation classified with PLS-DA.

Packages

windowsFonts(A = windowsFont("DejaVu Sans")) 
palette("Tableau10")

library(caret) #classification using train()
#library(kernlab)    
library(MASS) #generate the samples and calculate the augmented noise
library(dplyr) #data handling 
library(knitr) #knit options for the Rmarkdown report
library(kableExtra) #table options

Local sources files containing the functions

source("predict_with_uncertainty_parallel_v4.R")
source("plot_functions.R")

2D simulations with linear separation pattern

set.seed(1234)
n_points <- 10  # Number of points per class
n_replicates <- 5  # Number of replicates 
noise_sd <- 0.1  # Standard deviation 
rho <- 0 # Correlation coefficient

# Define the covariance matrix
cov_matrix <- matrix(c(noise_sd^2, rho * noise_sd^2, 
                       rho * noise_sd^2, noise_sd^2), nrow = 2)

# Function to generate 2D points with correlated noise
generate_class <- function(n_points, n_replicates, center_x, center_y, cov_matrix) {
  points <- matrix(0, nrow = n_points * n_replicates, ncol = 2)
  
  for (i in 1:n_points) {
    base_point <- c(center_x + rnorm(1, 0, 0.2), center_y + rnorm(1, 0, 0.2)) 
    
    for (j in 1:n_replicates) {
      noise <- mvrnorm(1, mu = c(0, 0), Sigma = cov_matrix)  # Correlated noise
      points[(i - 1) * n_replicates + j, ] <- base_point + noise 
    }
  }
  return(points)
}

class1 <- generate_class(n_points, n_replicates, center_x = 0.8, center_y = 1.6, cov_matrix)
class2 <- generate_class(n_points, n_replicates, center_x = 1.3, center_y = 1.2, cov_matrix)

data <- data.frame(
  Class = factor(rep(c("Class 1", "Class 2"), each = n_points * n_replicates)),
  Sample_ID = rep(1:(2 * n_points), each = n_replicates),
  X1 = c(class1[, 1], class2[, 1]),
  X2 = c(class1[, 2], class2[, 2])
)

# Mean center the data
data$X1 <- data$X1 - mean(data$X1)
data$X2 <- data$X2 - mean(data$X2)

# Calculate the mean and standard deviation of each sample
mean_data <- aggregate(cbind(X1, X2) ~ Class + Sample_ID, data = data, FUN = mean)
sd_data <- aggregate(cbind(X1, X2) ~ Class + Sample_ID, data = data, FUN = sd)

rm(class1, class2)

Fitting the PLS-DA model

pls_model <- train(y = mean_data$Class,
                   x = mean_data[,3:4], 
                   method = "pls")
pls_model
## Partial Least Squares 
## 
## 20 samples
##  2 predictor
##  2 classes: 'Class 1', 'Class 2' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 20, 20, 20, 20, 20, 20, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9473016  0.8804084
## 
## Tuning parameter 'ncomp' was held constant at a value of 1

Local uncertainty

The predict_with_uncertainty_parallel() function is used to estimate the uncertainty of each sample.

loc_pred <- predict_with_uncertainty(data[,-1],
                                                pls_model,
                                                n_sim = 500,
                                                return_simulations = T)

This function returns the mean and standard deviation of the \(\hat{y}\) of each sample based on the Monte Carlo simulation. It also returns the area above and bellow the decision boundary (in this case, zero) of each sample. The classification output indicates if the sample is being classified with certainty or not. For samples where the distribution crosses the decision boundary (considering the 95% confidence level), the sample will be assigned as uncertain. Once return_simulations = T, the augmented noise will be saved. The results for this 2D simulated data is displayed in Table 1.

Table 1. Results of the predict_with_uncertainty_parallel() function.
class sample y_mean y_sd area_above_0 area_below_0 classification
Class 1 1 0.5082 0.0466 1.000 0.000 Class 1
Class 1 2 0.4956 0.0573 1.000 0.000 Class 1
Class 1 3 0.1128 0.0484 0.990 0.010 Class 1
Class 1 4 0.4369 0.0425 1.000 0.000 Class 1
Class 1 5 0.2655 0.0605 1.000 0.000 Class 1
Class 1 6 0.6282 0.0535 1.000 0.000 Class 1
Class 1 7 0.4392 0.0579 1.000 0.000 Class 1
Class 1 8 0.2487 0.0610 1.000 0.000 Class 1
Class 1 9 0.7196 0.0728 1.000 0.000 Class 1
Class 1 10 0.2058 0.0357 1.000 0.000 Class 1
Class 2 1 -0.6100 0.0773 0.000 1.000 Class 2
Class 2 2 -0.2612 0.0701 0.000 1.000 Class 2
Class 2 3 -0.5392 0.0342 0.000 1.000 Class 2
Class 2 4 -0.7054 0.0492 0.000 1.000 Class 2
Class 2 5 -0.2287 0.0783 0.002 0.998 Class 2
Class 2 6 0.0412 0.1038 0.654 0.346 Uncertain
Class 2 7 -0.2905 0.0761 0.000 1.000 Class 2
Class 2 8 -0.5046 0.0653 0.000 1.000 Class 2
Class 2 9 -0.4913 0.0736 0.000 1.000 Class 2
Class 2 10 -0.4862 0.0590 0.000 1.000 Class 2

These results can be visualized, as showed in Figure 1.

Figure 1. Uncertainty of the classification outputs for the PLS-DA model. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.

Figure 1. Uncertainty of the classification outputs for the PLS-DA model. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.

Correlated noise

The previous examples showed how the uncertainty was estimated for Independent and Identically Distributed noise (iid noise). The next examples show how the uncertainty can be estimated for correlated noise, easily found in spectroscopy data.

Positive covariance

For this example, the covariance between the X1 and X2 is 0.8. The results are showed in Table 2.

. .
Table 2. Results of the predict_with_uncertainty_parallel() function for variables with positive covariance.
class sample y_mean y_sd area_above_0 area_below_0 classification
Class 1 1 0.5406 0.0315 1.000 0.000 Class 1
Class 1 2 0.4871 0.0401 1.000 0.000 Class 1
Class 1 3 0.1887 0.0144 1.000 0.000 Class 1
Class 1 4 0.5337 0.0211 1.000 0.000 Class 1
Class 1 5 0.2957 0.0349 1.000 0.000 Class 1
Class 1 6 0.5424 0.0346 1.000 0.000 Class 1
Class 1 7 0.4496 0.0193 1.000 0.000 Class 1
Class 1 8 0.1962 0.0177 1.000 0.000 Class 1
Class 1 9 0.6809 0.0238 1.000 0.000 Class 1
Class 1 10 0.2320 0.0237 1.000 0.000 Class 1
Class 2 1 -0.6341 0.0177 0.000 1.000 Class 2
Class 2 2 -0.2022 0.0341 0.000 1.000 Class 2
Class 2 3 -0.5983 0.0300 0.000 1.000 Class 2
Class 2 4 -0.6384 0.0142 0.000 1.000 Class 2
Class 2 5 -0.2666 0.0571 0.000 1.000 Class 2
Class 2 6 0.0045 0.0350 0.551 0.449 Uncertain
Class 2 7 -0.3333 0.0299 0.000 1.000 Class 2
Class 2 8 -0.5318 0.0307 0.000 1.000 Class 2
Class 2 9 -0.6027 0.0196 0.000 1.000 Class 2
Class 2 10 -0.3438 0.0281 0.000 1.000 Class 2

The results can also be visualized, as showed in Figure 2.

Figure 3. Uncertainty of the classification outputs for the PLS-DA model with positive covariance. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.

Figure 3. Uncertainty of the classification outputs for the PLS-DA model with positive covariance. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.

Negative Covariance

For this example, the covariance between the X1 and X2 is -0.8. The results are showed in Table 3 and can be visualized in Figure 3.

. .
Table 3. Results of the predict_with_uncertainty_parallel() function for variables with negative covariance.
class sample y_mean y_sd area_above_0 area_below_0 classification
Class 1 1 0.5719 0.0618 1.000 0.000 Class 1
Class 1 2 0.4288 0.0579 1.000 0.000 Class 1
Class 1 3 0.1413 0.0890 0.944 0.056 Uncertain
Class 1 4 0.4316 0.0902 1.000 0.000 Class 1
Class 1 5 0.2945 0.1363 0.985 0.015 Class 1
Class 1 6 0.5971 0.0729 1.000 0.000 Class 1
Class 1 7 0.4407 0.1167 1.000 0.000 Class 1
Class 1 8 0.2700 0.0712 1.000 0.000 Class 1
Class 1 9 0.6664 0.0696 1.000 0.000 Class 1
Class 1 10 0.1202 0.0653 0.967 0.033 Uncertain
Class 2 1 -0.5800 0.1028 0.000 1.000 Class 2
Class 2 2 -0.2759 0.0879 0.001 0.999 Class 2
Class 2 3 -0.5170 0.0757 0.000 1.000 Class 2
Class 2 4 -0.7555 0.0861 0.000 1.000 Class 2
Class 2 5 -0.2188 0.0956 0.011 0.989 Class 2
Class 2 6 0.0764 0.0591 0.902 0.098 Uncertain
Class 2 7 -0.2266 0.0691 0.001 0.999 Class 2
Class 2 8 -0.4275 0.0789 0.000 1.000 Class 2
Class 2 9 -0.4655 0.1157 0.000 1.000 Class 2
Class 2 10 -0.5930 0.1204 0.000 1.000 Class 2
Figure 3. Uncertainty of the classification outputs for the PLS-DA model with negative covariance. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.

Figure 3. Uncertainty of the classification outputs for the PLS-DA model with negative covariance. (a) 2D data simulation. Circles represents the replicates, while the triangles represents the mean of replicates. (b) Augmented noise generated through eigen decomposition of the original covariance matrix for each sample. 50 points are showed for each sample. (c) Estimated uncertainty for each sample.