Overlay graphics in ggplot2 considering fitted lines

Question

I am trying to superimpose two graphs made in ggplot2, in which the first contains the adjusted line using a polynomial model of degree 2 and the other contains the adjusted lines of a model using the gamm function, as can be seen in the images below.

What I am trying to accomplish is to overlay each of the information in both graphs on a single face, that is, an example would be as in the following assembly:

However, when trying to accomplish what I described earlier, I performed the following computational routine:

################################################################################
################################################################################
################################################################################
library(ggplot2)
library(splines)
library(nlme)
library(lattice)
library(latticeExtra)
require(gamm4)
library(RColorBrewer)
library(ggpmisc)
library(dplyr)
################################################################################
################################################################################
################################################################################
dados1 = read.table("CapilaridadeSemTempo0.csv", header = T, sep=";", dec=",")
dados2 <- reshape(cbind(id=1:nrow(dados1), dados1), 
                  varying=5:45, 
                  v.names="massaseca",
                  timevar="Tempo", 
                  times=as.numeric(gsub("X", "", tail(names(dados1), -3))), 
                  direction="long", sep="")
dados=dados2[order(dados2$id), ]
dados = na.omit(dados)
dados$Teor  <- factor(dados$Teor)
dados$Fator  <- factor(dados$Fator)
################################################################################
################################################################################
################################################################################
dados[dados$Teor == 0 & dados$Fator == "Am", "Trat"] = "Am1:0%"
dados[dados$Teor == 6 & dados$Fator == "Am", "Trat"] = "Am2:6%"
dados[dados$Teor == 8 & dados$Fator == "Am", "Trat"] = "Am3:8%"
dados[dados$Teor == 10 & dados$Fator == "Am", "Trat"] = "Am4:10%"
dados[dados$Teor == 12 & dados$Fator == "Am", "Trat"] = "Am5:12%"
dados[dados$Teor == 0 & dados$Fator == "As", "Trat"] = "As1:0%"
dados[dados$Teor == 6 & dados$Fator == "As", "Trat"] = "As2:6%"
dados[dados$Teor == 8 & dados$Fator == "As", "Trat"] = "As3:8%"
dados[dados$Teor == 10 & dados$Fator == "As", "Trat"] = "As4:10%"
dados[dados$Teor == 12 & dados$Fator == "As", "Trat"] = "As5:12%"
################################################################################
################################################################################
################################################################################
dados$Trat <- factor(dados$Trat)
################################################################################
################################################################################
################################## Model #######################################
################################################################################
################################################################################
fit4.gamm <- gamm(massaseca~factor(Trat)+s(Tempo,k=10,bs="ps",m=2,
                                           by=factor(Trat)),
                  random=list(id=pdSymm(~Tempo)),data=dados)
################################################################################
################################################################################
################################## Ajuste quadrático ###########################
################################################################################
x11()
my.formula <- y ~ x
ylim_sup <- 1.1 * max(dados$massaseca)
ylim_inf <- min(dados$massaseca)
shape_brks <- unique(dados$Trat)
shape_vals <- rep(1, 10)

label_y_npc <- rep(0.9, 10)
label_x_npc <- rep(0.91, 10)

p1 = dados %>%
  group_by(Tempo, Fator, Trat) %>%
  summarise(massaseca = mean(massaseca, na.rm = TRUE),.groups = 'drop') %>%
  ggplot(aes(x = Tempo, y = massaseca, shape = Trat,color = Trat)) +
  geom_point() + 
  stat_smooth(method = "lm", se = FALSE,
              formula = y ~ poly(x, 2, raw = TRUE),
              linetype = 1, 
              size = 1.1) + scale_shape_manual(name = "Trat", 
                                               breaks = shape_brks, 
                                               values = shape_vals) +
  coord_cartesian(xlim=c(0,1420)) + 
  stat_poly_eq(formula = y ~ poly(x, 2, raw = TRUE), 
               eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., sep = "*plain(\", \")~")),
               label.x.npc = label_x_npc,
               y = 0.8,  angle=90, 
               label.y.npc = label_y_npc,
               parse = TRUE, size = 4) +
  ylim(ylim_inf, ylim_sup) +
  labs(title = "",
       x = "Time (Minutes)",
       y = "Weight (mg)",
       color = "Trat") + 
  theme(legend.position = "none",axis.title = element_text(size = 23,color="black"),
        axis.text = element_text(size = 18,color="black"),
        text = element_text(size = 20,color="black")) + 
  facet_wrap(~Trat,ncol=5,nrow=2)
#################################################################################
################################# ggplot 2 - ajuste gamm ########################
#################################################################################
x11()

p2 = ggplot(dados, aes(x = Tempo, y = fitted(fit4.gamm$lme), group=id)) + 
  facet_wrap(~Trat,ncol=5,nrow=2) +
  xlab("Time (Minutes)") +
  geom_point(size=3, color="#969696") +
  geom_line(aes(x=Tempo,y=fitted(fit4.gamm$lme))) +
  ylab("Weight (mg)") +
  theme(legend.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 22,color="black"),
        axis.text.x = element_text(size = 18, colour = "black"),
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 25)) 

p1+p2

where I came across the following error:

Erro: Can't add `p2` to a ggplot object.
Run `rlang::last_error()` to see where the error occurred.

Please adjust your question to be minimal and reproducible, without external links. https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example — Jon Spring, Feb 20 '21 at 01:19
Can you just add `+geom_line(data = dados, aes(x=Tempo,y=fitted(fit4.gamm$lme, group = id))) +` to the first plot? that's a guess, can't confirm or help any further w/o reproducible example. — Jon Spring, Feb 20 '21 at 01:21
Hellow @JonSpring, thanks for the feedback, but, this link is the data. As for your comment, it was not useful since my interest is to superimpose both graphs in one as in the example. I thank! Note: your comment is already being used in the second part of the graph. — user55546, Feb 20 '21 at 13:19
This site discourages external links because 1) it creates a dependency on another site that might not be maintained into the future, reducing the future value of your question and potential answers, and 2) many potential answerers are untrusting of external links from strangers. — Jon Spring, Feb 20 '21 at 17:17
You can't add separate ggplot2 objects together like `p1+p2`. What you can do is add layers, inside the creation of `p1` like I suggested above. — Jon Spring, Feb 20 '21 at 17:20

user55546 · Accepted Answer · 2021-02-24T13:05:53.717

It is known that ggplot2 has its structure defined in layers, so the solution was structured to organize the syntax step by step in such a way that the desired specificities were incorporated.

Solution:

my.formula <- y ~ x
ylim_sup <- 1.1 * max(dados$massaseca)
ylim_inf <- min(dados$massaseca)
shape_brks <- unique(dados$Trat)
shape_vals <- rep(1, 10)

label_y_npc <- rep(1.0, 10)
label_x_npc <- rep(0.91, 10)

dados %>%
  group_by(Tempo, Fator, Trat) %>%
  summarise(massaseca = mean(massaseca, na.rm = TRUE),.groups = 'drop') %>%
  ggplot(aes(x = Tempo, y = massaseca,color = Trat)) +
  geom_point(data= dados, aes(x=Tempo,y=fitted(fit4.gamm$lme), group=id), color="#969696") +
  geom_line(data= dados, aes(x=Tempo,y=fitted(fit4.gamm$lme), group=id),
            col = "black") +
  stat_smooth(method = "lm", se = FALSE,
              formula = y ~ poly(x, 2, raw = TRUE),
              linetype = 1, 
              size = 0.7) + scale_shape_manual(name = "Trat", 
                                               breaks = shape_brks, 
                                               values = shape_vals) +
  geom_point() + 
  coord_cartesian(xlim=c(0,1420)) + 
  stat_poly_eq(formula = y ~ poly(x, 2, raw = TRUE), 
               eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., sep = "*plain(\", \")~")),
               label.x.npc = label_x_npc,
               y = 1.0,  angle=90, 
               label.y.npc = label_y_npc,
               parse = TRUE, size = 4.3) +
  ylim(ylim_inf, ylim_sup) +
  labs(title = "",
       x = "Time (Minutes)",
       y = "Weight (mg)",
       color = "Trat") + 
  theme(legend.position="none",legend.title = element_text(size = 20),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 22,color="black"),
        axis.text.x = element_text(size = 18, colour = "black"),
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title = element_text(size = 22),
        axis.text = element_text(size = 25),
        text = element_text(size = 20,color="black"))+ 
  facet_wrap(~Trat,ncol=5,nrow=2)

Overlay graphics in ggplot2 considering fitted lines

1 Answers1