Overlay graphics in ggplot2 considering adjusted lines

Asked

Viewed 92 times

1

I am trying to superimpose two graphs made in ggplot2, in which the first contains the line adjusted by means of 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’m trying to accomplish, is to overlay each of the information of both graphs on a single face, that is, an example would be:

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

Dice: https://drive.google.com/file/d/1QaBqcA7Y2F_eh6BOWP6Rlywzv6B4Mutl/view?usp=sharing

################################################################################
################################################################################
################################################################################
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

in which I came across the following mistake:

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

    Hello Breno. To answer the question use the answer field, I think you should know the system of the site https://answall.com/tour, I also recommend that you read: https://answall.com/help/self-answer. Thank you for understanding.

  • 1

    Provide a example minimum reproducible. This is not only important for those who try to help you, but also to help you locate and understand the problem in your code.

  • Hello @Guilhermenascimento, where is the answer field? Sorry for the question.

  • @Carloseduardolagosta I managed to accomplish the solution, I intend to post but how do I post as response?

  • 1

    Hello Breno. Scroll and you will see a text field with the title Sua resposta, write the report in a way that is useful to other users, explaining clearly what you did and finally click the blue button below written Publique sua resposta

1 answer

0


It is known that ggplot2 has its layered structure, so the solution was structured in organizing the syntax step by step in such a way that the desired specifics 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)

Browser other questions tagged

You are not signed in. Login or sign up in order to post.