考虑拟合线在ggplot2中叠加图形
Overlay graphics in ggplot2 considering fitted lines
我正在尝试叠加在 ggplot2 中制作的两个图形，其中第一个包含使用 2 次多项式模型调整的线，另一个包含使用 gamm 函数的模型的调整线，如图所示在下面的图片中。
我想要完成的是将两个图中的每个信息都叠加在一个面上，也就是说，一个示例如下所示：
但是，在尝试完成我之前描述的内容时，我执行了以下计算例程：
################################################################################
################################################################################
################################################################################
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
我遇到以下错误的地方：
Erro: Can't add `p2` to a ggplot object.
Run `rlang::last_error()` to see where the error occurred.
众所周知，ggplot2 的结构是按层定义的，因此解决方案的结构是逐步组织语法，以便合并所需的特性。
解决方案：
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)
我正在尝试叠加在 ggplot2 中制作的两个图形，其中第一个包含使用 2 次多项式模型调整的线，另一个包含使用 gamm 函数的模型的调整线，如图所示在下面的图片中。
我想要完成的是将两个图中的每个信息都叠加在一个面上，也就是说，一个示例如下所示：
但是，在尝试完成我之前描述的内容时，我执行了以下计算例程：
################################################################################
################################################################################
################################################################################
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
我遇到以下错误的地方：
Erro: Can't add `p2` to a ggplot object.
Run `rlang::last_error()` to see where the error occurred.
众所周知，ggplot2 的结构是按层定义的，因此解决方案的结构是逐步组织语法，以便合并所需的特性。
解决方案：
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)