With the following data, I am trying to create a plot as shown, but the order of legends as: mgCST Group, CST Group, Treatment Status. However, ggplot is ignoring all attempts to order the legends

data <- structure(list(ID = c(138L, 139L, 140L, 71L, 72L, 73L), UID = c("UAB005_W10D5", 
"UAB005_W10D6", "UAB005_W10D7", "UAB005_W1D1", "UAB005_W1D2", 
"UAB005_W1D3"), SID = c("UAB005 - CST", "UAB005 - CST", "UAB005 - CST", 
"UAB005 - mgCST", "UAB005 - mgCST", "UAB005 - mgCST"), SERIAL = c(68, 
69, 70, 1, 2, 3), WEEK = c(10, 10, 10, 1, 1, 1), DAY = c(5, 6, 
7, 1, 2, 3), FULL.STATUS = structure(c(5L, 5L, 5L, 2L, 2L, 2L
), levels = c("Prior to SBV", "Prior to MET", "BV DX", "During MET", 
"After MET"), class = "factor"), NUGENT_SCORE = c(7, 7, 7, 0, 
3, 5), NUGENT_CLASS = c("BV", "BV", "BV", "NO_BV", "NO_BV", "INTER_BV"
), mgCST_CST = structure(c(28L, 28L, 28L, NA, NA, NA), levels = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", 
"I", "III", "IV-A", "IV-B", "IV-C"), class = "factor"), subCST = c("IV-B", 
"IV-B", "IV-B", "III-B", "III-B", "III-B")), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -6L), groups = structure(list(
    SID = c("UAB005 - CST", "UAB005 - mgCST"), .rows = structure(list(
        1:3, 4:6), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -2L), .drop = TRUE))

MET <- ggplot(mapping = aes(x = SERIAL, y = SID,  shape = FULL.STATUS)) +
              geom_point(data = data.CST, aes(color = CST), size = 2) +
              scale_color_manual(values = CST.colors,
                                 name = "CST Group") +
              ggnewscale::new_scale_color() +
              geom_point(data = data.mgCST, aes(color = mgCST), size = 2) +
              scale_color_manual(values = mgCST.colors,
                                 name = "mgCST Group") +
              scale_shape_manual(values = status.shapes,
                                 name = "Treatment Status") +
              guides(color = guide_legend(order = 1)) +
              xlim(c(0, 70)) +
              xlab("Study Day") +
              ylab("SID") +
              theme_bw() + 
              theme(text = element_text(size = 12)) + 
              ggtitle("BV Dynamics MET Study")

I have tried all iterations of the guide function:

MET <- ggplot(mapping = aes(x = SERIAL, y = SID,  shape = FULL.STATUS)) +
              geom_point(data = data.CST, aes(color = CST), size = 2) +
              scale_color_manual(values = CST.colors,
                                 name = "CST Group") +
              ggnewscale::new_scale_color() +
              geom_point(data = data.mgCST, aes(color = mgCST), size = 2) +
              scale_color_manual(values = mgCST.colors,
                                 name = "mgCST Group") +
              scale_shape_manual(values = status.shapes,
                                 name = "Treatment Status") +
              guides(color = guide_legend(order = 1),
                     color = guide_legend(order = 2),
                     shape = guide_legend(order = 3)) +
              xlim(c(0, 70)) +
              xlab("Study Day") +
              ylab("SID") +
              theme_bw() + 
              theme(text = element_text(size = 12)) + 
              ggtitle("BV Dynamics MET Study")

MET <- ggplot(mapping = aes(x = SERIAL, y = SID,  shape = FULL.STATUS)) +
              geom_point(data = data.CST, aes(color = CST), size = 2) +
              scale_color_manual(values = CST.colors,
                                 name = "CST Group") +
              ggnewscale::new_scale_color() +
              geom_point(data = data.mgCST, aes(color = mgCST), size = 2) +
              scale_color_manual(values = mgCST.colors,
                                 name = "mgCST Group") +
              scale_shape_manual(values = status.shapes,
                                 name = "Treatment Status") +
              xlim(c(0, 70)) +
              xlab("Study Day") +
              ylab("SID") +
              theme_bw() + 
              theme(text = element_text(size = 12)) + 
              ggtitle("BV Dynamics MET Study")

MET + guides(color = guide_legend(order = 1),
             shape = guide_legend(order = 3))

as well as the following:

guides_merge <- function(gdefs) {
  gdefs <- lapply(gdefs, function(g) { g$hash <- paste(g$order, g$hash, sep = "z"); g})
  tapply(gdefs, sapply(gdefs, function(g)g$hash), function(gs)Reduce(guide_merge, gs))
}
environment(guides_merge) <- environment(ggplot)
assignInNamespace("guides_merge", guides_merge, pos = "package:ggplot2")

MET <- ggplot(mapping = aes(x = SERIAL, y = SID,  shape = FULL.STATUS)) +
              geom_point(data = data.CST, aes(color = CST), size = 2) +
              scale_color_manual(values = CST.colors,
                                 name = "CST Group") +
              ggnewscale::new_scale_color() +
              geom_point(data = data.mgCST, aes(color = mgCST), size = 2) +
              scale_color_manual(values = mgCST.colors,
                                 name = "mgCST Group") +
              scale_shape_manual(values = status.shapes,
                                 name = "Treatment Status") +
              guides(color = guide_legend(order = 1),
                     color = guide_legend(order = 2),
                     shape = guide_legend(order = 3)) +
              xlim(c(0, 70)) +
              xlab("Study Day") +
              ylab("SID") +
              theme_bw() + 
              theme(text = element_text(size = 12)) + 
              ggtitle("BV Dynamics MET Study")

I have updated ggplot2 and would appreciate any help!