# Your data could look like this:
@@ -170,13 +170,19 @@ package.
# select AB results dynamically
select(mo, aminoglycosides(), betalactams()) %>%
# replace NAs with NI (not-interpretable)
- mutate(across(where(is.sir),
- ~replace_na(.x, "NI")),
- # make factors of SIR columns
- across(where(is.sir),
- as.integer),
- # get Gramstain of microorganisms
- mo = as.factor(mo_gramstain(mo))) %>%
+ mutate(
+ across(
+ where(is.sir),
+ ~ replace_na(.x, "NI")
+ ),
+ # make factors of SIR columns
+ across(
+ where(is.sir),
+ as.integer
+ ),
+ # get Gramstain of microorganisms
+ mo = as.factor(mo_gramstain(mo))
+ ) %>%
# drop NAs - the ones without a Gramstain (fungi, etc.)
drop_na()
#> ℹ For `aminoglycosides()` using columns GEN (gentamicin), TOB (tobramycin), AMK
@@ -334,7 +340,7 @@ performance.
set.seed(123) # For reproducibility
data_split <- initial_split(data, prop = 0.8) # 80% training, 20% testing
training_data <- training(data_split) # Training set
-testing_data <- testing(data_split) # Testing set
+testing_data <- testing(data_split) # Testing set
# Fit the workflow to the training data
fitted_workflow <- resistance_workflow %>%
@@ -354,7 +360,7 @@ called since they are stored in the recipe.
# Make predictions on the testing set
predictions <- fitted_workflow %>%
- predict(testing_data) # Generate predictions
+ predict(testing_data) # Generate predictions
probabilities <- fitted_workflow %>%
predict(testing_data, type = "prob") # Generate probabilities
@@ -529,8 +535,8 @@ predictors can easily and agnostically selected using the new
# Define the recipe
mic_recipe <- recipe(esbl ~ ., data = training_data) %>%
- remove_role(genus, old_role = "predictor") %>% # Remove non-informative variable
- step_mic_log2(all_mic_predictors()) # Log2 transform all MIC predictors
+ remove_role(genus, old_role = "predictor") %>% # Remove non-informative variable
+ step_mic_log2(all_mic_predictors()) # Log2 transform all MIC predictors
prep(mic_recipe)
#>
@@ -665,27 +671,38 @@ status.
ggplot(predictions, aes(x = esbl, fill = .pred_class)) +
geom_bar(position = "stack") +
- labs(title = "Predicted vs Actual ESBL Status",
- x = "Actual ESBL",
- y = "Count") +
+ labs(
+ title = "Predicted vs Actual ESBL Status",
+ x = "Actual ESBL",
+ y = "Count"
+ ) +
theme_minimal()
And plot the certainties too - how certain were the actual
predictions?
predictions %>%
- mutate(certainty = ifelse(.pred_class == "FALSE",
- .pred_FALSE,
- .pred_TRUE),
- correct = ifelse(esbl == .pred_class, "Right", "Wrong")) %>%
- ggplot(aes(x = seq_len(nrow(predictions)),
- y = certainty,
- colour = correct)) +
- scale_colour_manual(values = c(Right = "green3", Wrong = "red2"),
- name = "Correct?") +
+ mutate(
+ certainty = ifelse(.pred_class == "FALSE",
+ .pred_FALSE,
+ .pred_TRUE
+ ),
+ correct = ifelse(esbl == .pred_class, "Right", "Wrong")
+ ) %>%
+ ggplot(aes(
+ x = seq_len(nrow(predictions)),
+ y = certainty,
+ colour = correct
+ )) +
+ scale_colour_manual(
+ values = c(Right = "green3", Wrong = "red2"),
+ name = "Correct?"
+ ) +
geom_point() +
- scale_y_continuous(labels = function(x) paste0(x * 100, "%"),
- limits = c(0.5, 1)) +
+ scale_y_continuous(
+ labels = function(x) paste0(x * 100, "%"),
+ limits = c(0.5, 1)
+ ) +
theme_minimal()
# Define the recipe
resistance_recipe_time <- recipe(res_AMX ~ year + gramstain, data = data_time) %>%
- step_dummy(gramstain, one_hot = TRUE) %>% # Convert categorical to numerical
- step_normalize(year) %>% # Normalise year for better model performance
- step_nzv(all_predictors()) # Remove near-zero variance predictors
+ step_dummy(gramstain, one_hot = TRUE) %>% # Convert categorical to numerical
+ step_normalize(year) %>% # Normalise year for better model performance
+ step_nzv(all_predictors()) # Remove near-zero variance predictors
resistance_recipe_time
#>
@@ -929,9 +951,11 @@ sets.
ggplot(predictions_time, aes(x = year)) +
geom_point(aes(y = res_AMX, color = "Actual")) +
geom_line(aes(y = .pred, color = "Predicted")) +
- labs(title = "Predicted vs Actual AMX Resistance Over Time",
- x = "Year",
- y = "Resistance Proportion") +
+ labs(
+ title = "Predicted vs Actual AMX Resistance Over Time",
+ x = "Year",
+ y = "Resistance Proportion"
+ ) +
theme_minimal()
Additionally, we can visualise resistance trends in
@@ -939,13 +963,17 @@ sets.
ggplot(data_time, aes(x = year, y = res_AMX, color = gramstain)) +
geom_line() +
- labs(title = "AMX Resistance Trends",
- x = "Year",
- y = "Resistance Proportion") +
+ labs(
+ title = "AMX Resistance Trends",
+ x = "Year",
+ y = "Resistance Proportion"
+ ) +
# add a linear model directly in ggplot2:
- geom_smooth(method = "lm",
- formula = y ~ x,
- alpha = 0.25) +
+ geom_smooth(
+ method = "lm",
+ formula = y ~ x,
+ alpha = 0.25
+ ) +
theme_minimal()