How to build a factors mean graphs (multifactor ANOVA) in ggplot2

Question

I am performing a L9 orthogonal design, and to see visually how the different factor levels are, I am interested in building a means plot. With the code:

plot.design(PCR ~ Temperature+Time+Ratio+Pretreatment, data = scenedesmus)

I obtain:

I would like to improve this graph with ggplot2, but I was not capable of it. Moreover, it would be great if the different levels for each factor were slightly separated from each other, with different colors per factor, and with the standard deviation of each mean factor incorporated. Alone I obtain each of those graphs in ggplot2 individually, but I didn´t know how to merge them:

ggline(scenedesmus, x = "Temperature", y = "PCR",
       add = "mean_sd", size=1,
       ylab ="Protein-carbohydrate ratio", legend="right", error.plot = "errorbar") + theme(axis.text=element_text(size=20), axis.title=element_text(size=24,face="bold"), legend.title=element_text(size=24, face="bold"),legend.text=element_text(size=22))+     scale_x_discrete(labels = scales::label_parse())

I would like to obtain something like this:

My dataset is the following:

> scenedesmus
Temperature Time Ratio Pretreatment       PRY      CRY
1           20  0.5     3         None  7.106190 12.99137
2           20  0.5     3         None  6.991073 13.20371
3           20    1     6       Mortar  9.816545 14.62239
4           20    1     6       Mortar 10.093768 14.41567
5           20    2    12        Discs 15.887290 20.85106
6           20    2    12        Discs 16.514740 21.13347
7           30  0.5     6        Discs 15.608507 20.75174
8           30  0.5     6        Discs 15.890457 20.37846
9           30    1    12         None  9.851556 13.19030
10          30    1    12         None 10.329157 12.74816
11          30    2     3       Mortar  9.815574 14.37999
12          30    2     3       Mortar 10.177421 15.15487
13          40  0.5    12       Mortar 12.097258 16.36536
14          40  0.5    12       Mortar 11.135055 17.34924
15          40    1     3        Discs 14.759191 22.44141
16          40    1     3        Discs 14.884651 22.43402
17          40    2     6         None  9.476980 14.08952
18          40    2     6         None 10.832856 16.03889
         PCR
1  0.5469931
2  0.5294781
3  0.6713366
4  0.7001942
5  0.7619416
6  0.7814496
7  0.7521542
8  0.7797673
9  0.7468787
10 0.8102468
11 0.6825854
12 0.6715611
13 0.7391992
14 0.6418181
15 0.6576767
16 0.6634856
17 0.6726260
18 0.6754117

Answer 1

You would need to pivot all your factors into a single text column and plot them on the x axis.

library(tidyverse)

scenedesmus %>%
  mutate(across(Temperature:Pretreatment, as.character)) %>%
  pivot_longer(Temperature:Pretreatment) %>%
  mutate(value = factor(value, levels(factor(value))[
    c(10:12, 6, 9, 3, 5, 7:8, 1:2, 4)])) %>%
  ggplot(aes(value, PCR, group = name, color = name)) +
  geom_point(stat = 'summary', fun = mean) +
  geom_errorbar(stat = 'summary', width = 0.1, alpha = 0.5) +
  geom_line(stat = 'summary', fun = mean) +
  geom_text(stat = 'summary', fun = mean, color = 'black',
            aes(label = after_stat(round(y, 2))),
            position = position_nudge(0.4, 0.01)) +
  facet_grid(~name, scales = 'free_x', switch = 'x') +
  scale_color_brewer(palette = 'Set1', guide = 'none') +
  coord_cartesian(clip = 'off') +
  geom_vline(data = data.frame(a = 0.4, name = 'Pretreatment'),
             aes(xintercept = a)) +
  theme_minimal(base_size = 20) +
  theme(strip.placement = 'outside',
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.spacing.x = unit(0, 'mm'),
        axis.ticks = element_line(),
        axis.line.x = element_line())