R: Solubility Regression with Linear and Random Forest Models

Table 1: Well-known compounds in the dataset

Dataset Description

The original report published a dataset of compounds represented as SMILES strings. SMILES strings are a compact and text-based method of specifying chemical structures. This study will use a preprocessed dataset mentioned on moleculenet.ai and distributed by deepchem.io, which contains features parsed from these SMILES strings. You can browse the file on GitHub.This study uses a subset of these preprocessed features, which are listed in Table 2.

Exploratory visualization

df <- as_tibble(read.csv("data/delaney-processed.csv")) %>%
  select(
    compound = Compound.ID, 
    mw = Molecular.Weight, 
    h_bond_donors = Number.of.H.Bond.Donors, 
    rings = Number.of.Rings, 
    rotatable_bonds = Number.of.Rotatable.Bonds, 
    psa = Polar.Surface.Area, 
    solubility = measured.log.solubility.in.mols.per.litre
)

Before I dive into the machine learning model, let’s examine exploratory plots to get a feel for the data distribution. Figure 1 has histograms (for continuous variables) and bar plots (for discrete variables) to demonstrate the dataset’s values’ distributions. Figure 1a, 1b, 1d, 1e, and 1f show distributions of values favoring their respective range’s low ends. Solubility, our response variable, has a broader spread above and below its mean of -3.05.

p1 <- ggplot(df, aes(x = mw)) +
  geom_histogram(bins = 10) +
  labs(title = "(a)") +
  theme_minimal()

p2 <- ggplot(df, aes(x = psa)) +
  geom_histogram(bins = 10) +
  labs(title = "(b)") +
  theme_minimal()

p3 <- ggplot(df, aes(x = solubility)) +
  geom_histogram(bins = 10) +
  labs(title = "(c)") +
  theme_minimal()

p4 <- ggplot(df, aes(x = h_bond_donors)) +
  geom_bar() +
  labs(title = "(d)") +
  theme_minimal()

p5 <- ggplot(df, aes(x = rings)) +
  geom_bar() +
  labs(title = "(e)") +
  theme_minimal()

p6 <- ggplot(df, aes(x = rotatable_bonds)) +
  geom_bar() +
  labs(title = "(f)") +
  theme_minimal()

grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Figure 1: Histograms and bar plots of variables

Figure 1a, 1b, 1d, 1e, and 1f show distributions of values favoring their respective range’s low ends. Solubility, in Figure 1c, has a broader spread above and below its mean of -3.05.

A number of the features of the molecules require lots of atoms. For example, a five-ring molecule will likely have a higher molecular weight than a three-ring molecule. Molecular weight has a special relationship with all other variables that denote each molecule’s increasing structural complexity. Figure 2 plots molecular weight against all other variables, with a trend line for each relationship, as shown in the jittered scatter plots below.

alpha = 0.1
p1 <- ggplot(df, aes(x = psa, y = mw)) +
  geom_jitter(alpha = alpha) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "(a)") +
  theme_minimal()

p2 <- ggplot(df, aes(x = solubility, y = mw)) +
  geom_jitter(alpha = alpha) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "(b)") +
  theme_minimal()

p3 <- ggplot(df, aes(x = h_bond_donors, y = mw)) +
  geom_jitter(alpha = alpha, width = 0.1) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "(c)") +
  theme_minimal()

p4 <- ggplot(df, aes(x = rings, y = mw)) +
  geom_jitter(alpha = alpha, width = 0.1) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "(d)") +
  theme_minimal()

p5 <- ggplot(df, aes(x = rotatable_bonds, y = mw)) +
  geom_jitter(alpha = alpha, width = 0.1) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "(e)") +
  theme_minimal()

grid.arrange(p1, p2, p3, p4, p5, nrow = 3)

Figure 2: Relationship of molecular weight to other variables

Figures 2a, 2c, 2d, 2e all exhibit increasing molecular weight with increased structural complexity. Figure 2b stands out: in general, as molecular weight increases, solubility decreases.

Train/test split

All models use the same randomized train/test split. First, I shuffled the rows of the original dataset. Then, I selected the first 846 rows for the training set and the last 282 rows for the test dataset. I seeded the random number generator with a constant to ensure the same shuffles between model runs.

solubilityTrainTestSplit <- function(all_data, split_fraction = 0.75, random_seed = 13) {
  # Shuffle based on random seed
  set.seed(random_seed)
  sample_indecies <- sample(nrow(all_data), nrow(all_data))
  shuffled <- all_data[sample_indecies, ]
  
  # Train test split
  train_row <- round(nrow(shuffled) * split_fraction)
  test_row <- train_row + 1
  train <- shuffled[1:train_row, ]
  test <- shuffled[test_row:nrow(shuffled), ]
  
  # Now create the final list that is returned
  list(
    train = train,
    test = test
  )
}

split <- solubilityTrainTestSplit(df)
test <- split$test
train <- split$train

Table 3 is the first few rows of the train dataset:

Formula

All models use the same formula based on all the predictors.

formula <- solubility ~ mw + h_bond_donors + rings + rotatable_bonds + psa

Linear Model

Our first stop is the linear model. Table 3 shows the performance of the linear model after training when I use it to predict the solubilties in the test dataset.

analyzeSolubilityLinearModel <- function(model, test_data) {
  # Run the prediction
  predicted_solbilities <- predict(model, test_data)
  
  # Assemble the compounds and predictions back onto the
  # train features
  
  test_results <- test_data %>%
    mutate(prediction = predicted_solbilities) %>%
    mutate(residual = solubility - prediction)
  
  # Calculate the standard deviation and RMSE
  sd_solubilities <- sd(test_results$solubility)
  rmse <- sqrt(mean(test_results$residual ^ 2))
  
  # calculate r_squared
  rss <- sum(test_results$residual ^ 2)
  total_error <- test_results$solubility - mean(test_results$solubility)
  tss <- sum(total_error ^ 2)
  r_squared <- 1 - (rss / tss)
  
  # Create a list to return to the caller
  list(
    sd_solubilities = sd_solubilities,
    rmse = rmse,
    test_results = test_results,
    r_squared = r_squared
  )
}

lm_model <- lm(solubility ~ mw + h_bond_donors + rings + rotatable_bonds + psa, train)
lm_result <- analyzeSolubilityLinearModel(lm_model, test)
lm_result_df <- tibble(
  metric = c("RMSE", "R^2"),
  value = c(lm_result$rmse, lm_result$r_squared)
)
knitr::kable(lm_result_df)

Random Forest Regression

The second model is the random forest model. For the random forest, I use the ranger package. The model uses 500 trees and is trained and tested with the same train and test data as the original linear model. Table 4 shows the random forest’s prediction performance on the test data after training.

analyzeSolubilityRandomForestModel <- function(model, test_data) {
  # Run the prediction
  predicted_solbilities <- predict(model, test_data)$predictions
  
  # Assemble the compounds and predictions back onto the
  # train features
  
  test_results <- test_data %>%
    mutate(prediction = predicted_solbilities) %>%
    mutate(residual = solubility - prediction)
  
  # Calculate the standard deviation and RMSE
  sd_solubilities <- sd(test_results$solubility)
  rmse <- sqrt(mean(test_results$residual ^ 2))
  
  # calculate r_squared
  rss <- sum(test_results$residual ^ 2)
  total_error <- test_results$solubility - mean(test_results$solubility)
  tss <- sum(total_error ^ 2)
  r_squared <- 1 - (rss / tss)
  
  # Create a list to return to the caller
  list(
    sd_solubilities = sd_solubilities,
    rmse = rmse,
    test_results = test_results,
    r_squared = r_squared
  )
}

rf_model <- ranger(
  formula,
  train,
  num.trees = 500,
  respect.unordered.factors = "order"
)

rf_result <- analyzeSolubilityRandomForestModel(rf_model, test)
rf_result_df = tibble(
  metric = c("RMSE", "R^2"),
  value = c(rf_result$rmse, rf_result$r_squared)
)

knitr::kable(rf_result_df)

Comparing linear and random forest models

Table 5 shows the performance metrics of the linear and random forest models.

Table 5

model_comparison_df <- tibble(
  model = c("Linear", "Random forest"),
  rmse = c(lm_result$rmse, rf_result$rmse),
  r_squared = c(lm_result$r_squared, rf_result$r_squared)
)

knitr::kable(model_comparison_df)

model	rmse	r_squared
Linear	1.1955559	0.6721196
Random forest	0.8634554	0.8289766

Examining Table 5 shows reveals two comparisons between the models that are immediately apparent. The first comparison is their respective RMSE values when evaluated on the training dataset. Lower RMSE values are better. The standard deviation of log solubility in the test dataset is 2.092. This value compares nicely with the random forest’s RMSE of 0.867, but it does not compare well with the linear model’s RMSE of 1.196. The clear winner here, unsurprisingly, the random forest. The second comparison is their respective R^2 values. R^2 values closer to 1.0 are better. The R^2 value of the linear model is a lowly 0.672, while the R^2 value for the random forest is much better at 0.828.

Figure 3, which plots actual versus predicted solubilities, is my favorite comparison between the linear and random forest models. Recall that all solubilities are log solubilities, hence the negative values on the axes.

alpha <- 0.5

linear_result_plot <- ggplot(lm_result$test_results, aes(x = prediction, y = solubility)) +
  geom_jitter(alpha = alpha, width = 0.1) +
  stat_smooth(method = "lm", se = FALSE) +
  labs(title = "(a) Linear") +
  theme_minimal()

rf_result_plot <- ggplot(rf_result$test_results, aes(x = prediction, y = solubility)) +
  geom_jitter(alpha = alpha, width = 0.1) +
  stat_smooth(method = "lm", se = FALSE) +
  labs(title = "(b) Random forest") +
  theme_minimal()

grid.arrange(linear_result_plot, rf_result_plot, nrow = 1)

Figure 3: actual vs. predicted values for both models

In Figure 3a, the linear model fit is lacking in the area of high predicted solubility. In this area, the actual solubilities take a wide range, from relatively insoluble to very soluble. At the edge of low predicted solubility, we find a range of actual solubilities.

In Figure 3b, the random forest model trends slightly more reliable. The random forest predicts high solubilities where there are actual high solubilities. Near the low end of the actual solubilities, there are numerous under and over predictions from the random forest.

Conclusion

The random forest performs better than the linear model but still has deficiencies around molecular features causing over or under predicted solubilities.