Removing inside lines in multiple factor R plot

I am representing four categorical factors at three levels for a response variable. However, when I try to set up a border it creates a border for each factor, and I don´t know how to remove those veritcal lines.

Here is a MWE:

library(ggplot2)
library("scales")
library(ggpubr)
library(car)
library(tidyverse)
scenedesmus$Milling<-as.character(scenedesmus$Milling)
scenedesmus$Temperature<-factor(scenedesmus$Temperature)
scenedesmus$Time<-factor(scenedesmus$Time)
scenedesmus$Ratio<-factor(scenedesmus$Ratio)


dput(scenedesmus)
structure(list(Temperature = c(20, 20, 20, 20, 20, 20, 30, 30, 
30, 30, 30, 30, 40, 40, 40, 40, 40, 40), Time = c(0.5, 0.5, 1, 
1, 2, 2, 0.5, 0.5, 1, 1, 2, 2, 0.5, 0.5, 1, 1, 2, 2), Ratio = c(3, 
3, 6, 6, 12, 12, 6, 6, 12, 12, 3, 3, 12, 12, 3, 3, 6, 6), Milling = structure(c(1L, 
1L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 1L, 
1L), levels = c("None", "Mortar", "Discs"), class = "factor"), 
    PRY = c(7.10618979550317, 6.99107348052751, 9.81654489395678, 
    10.0937678454159, 15.8872899104855, 16.5147395153748, 15.6085073784574, 
    15.8904572330355, 9.85155639002801, 10.3291566375677, 9.81557388225615, 
    10.1774212169006, 12.0972576247432, 11.1350551614397, 14.7591913822601, 
    14.8846506719242, 9.47697977090569, 10.8328555963545), CRY = c(12.9913707456184, 
    13.2037056981015, 14.6223886369729, 14.4156689100426, 20.8510599220091, 
    21.1334682925674, 20.7517385553227, 20.3784601114164, 13.1903022986714, 
    12.7481614338955, 14.3799945987187, 15.1548695641213, 16.3653561008515, 
    17.3492383422838, 22.4414097199122, 22.4340213280367, 14.0895227253865, 
    16.0388931794408), PCR = c(0.546993072143667, 0.529478135939726, 
    0.671336615218633, 0.700194205929921, 0.761941597689038, 
    0.781449560798454, 0.752154203217537, 0.779767320305684, 
    0.746878742196859, 0.810246770966047, 0.682585366418063, 
    0.671561122571146, 0.739199168670324, 0.641818098394623, 
    0.657676659642481, 0.663485625438106, 0.672626032522028, 
    0.675411668071984)), row.names = c(NA, -18L), class = c("tbl_df", 
"tbl", "data.frame"))


scenedesmus %>%
  mutate(across(Temperature:Milling, as.character)) %>%
  pivot_longer(Temperature:Milling) %>%
  mutate(value = factor(value, levels(factor(value))[
    c(12:10, 6, 9, 3, 5, 7:8, 1:2, 4)])) %>%
  ggplot(aes(value, PCR, group = name, color = name)) +
  geom_point(stat = 'summary', fun = mean, size=3) +
  geom_hline(yintercept = mean(scenedesmus$PCR, na.rm=TRUE),linetype='dotted', col = 'grey', size=1.5)+
  scale_y_continuous(limits = c(0.5, 0.9)) +
  geom_line(stat = 'summary', fun = mean, size=1) +
  facet_grid(~name, scales = 'free_x', switch = 'x') +
  scale_color_manual(values = c("#0072B2", "#D55E00", "#CC79A7","#009E73")) +
  coord_cartesian(clip = 'off') +
  geom_vline(data = data.frame(a = 0.4, name = 'Milling'),
             aes(xintercept = a)) +
  theme_classic(base_size = 20) +
  theme(strip.placement = 'outside',
        legend.position="none",
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.background = element_rect(color="black"),
        panel.spacing.x = unit(0, 'mm'),
        axis.ticks = element_line(),
        axis.line.x = element_line(),
        axis.title.y = element_text(size=20, face="bold"),
        strip.text = element_text(face = "bold")) +
  labs(y = "response")

Removing panel.background = element_rect(color="black") I obtained a plot but with the right and top borders missing

The issue is that panel.background will draw an outline around each panel. Instead, one option would be to fake your right and top border lines using geom_h/vline, i.e. add a geom_hline to each panel for which you set yintercept=Inf and a geom_vline for the rightmost panel.

library(scales)
library(car)
library(tidyverse)

scenedesmus %>%
  mutate(across(Temperature:Milling, as.character)) %>%
  pivot_longer(Temperature:Milling) %>%
  mutate(value = factor(value, levels(factor(value))[
    c(12:10, 6, 9, 3, 5, 7:8, 1:2, 4)
  ])) %>%
  ggplot(aes(value, PCR, group = name, color = name)) +
  geom_point(stat = "summary", fun = mean, size = 3) +
  geom_hline(yintercept = mean(scenedesmus$PCR, na.rm = TRUE), linetype = "dotted", col = "grey", size = 1.5) +
  geom_hline(
    yintercept = Inf, linewidth = 1
  ) +
  geom_vline(
    data = data.frame(x = Inf, name = "Time"),
    aes(xintercept = x), linewidth = 1
  ) +
  scale_y_continuous(limits = c(0.5, 0.9)) +
  geom_line(stat = "summary", fun = mean, size = 1) +
  facet_grid(~name, scales = "free_x", switch = "x") +
  scale_color_manual(values = c("#0072B2", "#D55E00", "#CC79A7", "#009E73")) +
  coord_cartesian(clip = "off") +
  geom_vline(
    data = data.frame(a = 0.4, name = "Milling"),
    aes(xintercept = a)
  ) +
  theme_classic(base_size = 20) +
  theme(
    strip.placement = "outside",
    legend.position = "none",
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.spacing.x = unit(0, "mm"),
    axis.ticks = element_line(),
    axis.line.x = element_line(),
    axis.title.y = element_text(size = 20, face = "bold"),
    strip.text = element_text(face = "bold")
  ) +
  labs(y = "response")