用误差条绘制线性回归分析的可变性

Question

我想制作类似于图 1 的图（来源：link）

在图 1 中，他们绘制了具有一年产量变化的回归分析。就我而言，我想绘制每个治疗组的两个位置和 4 个块之间的变异性。所以我想要的情节将具有三个方面的因素B.glucosidase、Protein、POX.C 的variable 和四种颜色的治疗因素。此外，在我目前的情节中，我有阻止和治疗的传说。我应该只接受治疗，因为该块应该用于制作可变性的误差线。

我尝试使用此代码，但显然不能满足我的要求。 （df.melted 的数据如下。）

ggplot(df.melted, aes(x = value, y = yield, color = as.factor(treatment))) + 
  geom_point(aes(shape= as.factor(block))) +
  stat_smooth(method = "lm", formula = y ~ x, col = "darkslategrey", se=F) +
  stat_poly_eq(formula = y~x, 
               # aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               aes(label =  ..rr.label..), 
               parse = TRUE) + 
  theme_classic() +
  geom_errorbar(aes(ymax = df.melted$yield+sd(df.melted$yield), ymin = df.melted$yield-sd(df.melted$yield)), width = 0.05)+
facet_wrap(~variable) 

数据：

df.melted <- structure(list(Location = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("M", "U"), class = "factor"), 
    treatment = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 
    3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
    2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 
    2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("CC", 
    "CCS", "CS", "SCS"), class = "factor"), block = c(1L, 2L, 
    3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
    2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 
    1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 
    4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 
    3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
    2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 
    1L, 2L, 3L, 4L), yield = c(5156L, 5157L, 5551L, 5156L, 4804L, 
    4720L, 4757L, 5021L, 4826L, 4807L, 4475L, 4596L, 4669L, 4588L, 
    4542L, 4592L, 5583L, 5442L, 5693L, 5739L, 5045L, 4902L, 5006L, 
    5086L, 4639L, 4781L, 4934L, 4857L, 4537L, 4890L, 4842L, 4608L, 
    5156L, 5157L, 5551L, 5156L, 4804L, 4720L, 4757L, 5021L, 4826L, 
    4807L, 4475L, 4596L, 4669L, 4588L, 4542L, 4592L, 5583L, 5442L, 
    5693L, 5739L, 5045L, 4902L, 5006L, 5086L, 4639L, 4781L, 4934L, 
    4857L, 4537L, 4890L, 4842L, 4608L, 5156L, 5157L, 5551L, 5156L, 
    4804L, 4720L, 4757L, 5021L, 4826L, 4807L, 4475L, 4596L, 4669L, 
    4588L, 4542L, 4592L, 5583L, 5442L, 5693L, 5739L, 5045L, 4902L, 
    5006L, 5086L, 4639L, 4781L, 4934L, 4857L, 4537L, 4890L, 4842L, 
    4608L), variable = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("B.glucosidase", 
    "Protein", "POX.C"), class = "factor"), value = c(1.600946, 
    1.474084, 1.433078, 1.532492, 1.198667, 1.193193, 1.214941, 
    1.360981, 1.853056, 1.690117, 1.544357, 1.825132, 1.695409, 
    1.764123, 1.903743, 1.538684, 0.845077, 1.011463, 0.857032, 
    0.989803, 0.859022, 0.919467, 1.01717, 0.861689, 0.972332, 
    0.952922, 0.804431, 0.742634, 1.195837, 1.267285, 1.08571, 
    1.20097, 6212.631579, 5641.403509, 4392.280702, 7120.701754, 
    5305.964912, 4936.842105, 5383.157895, 6077.894737, 5769.122807, 
    5016.842105, 5060.350877, 5967.017544, 5576.842105, 5174.035088, 
    5655.438596, 5468.77193, 7933.333333, 7000, 6352.982456, 
    8153.684211, 6077.894737, 4939.649123, 5002.807018, 6489.122807, 
    4694.035088, 5901.052632, 4303.859649, 6768.421053, 6159.298246, 
    6090.526316, 4939.649123, 5262.45614, 810.3024, 835.5242, 
    856.206, 759.8589, 726.2298, 792.6472, 724.7165, 699.3266, 
    500.9153, 634.8698, 637.9536, 648.8814, 641.0357, 623.3822, 
    555.2834, 520.8119, 683.3528, 595.9173, 635.4315, 672.4234, 
    847.2944, 745.5665, 778.3548, 735.8141, 395.2647, 570.4148, 
    458.0383, 535.3851, 678.0293, 670.7419, 335.2923, 562.5674
    )), row.names = c(NA, -96L), class = "data.frame")

Answer 1

library(dplyr)
library(ggplot2)
library(ggpmisc)

汇总数据框（这也可以使用stat_summary() 来完成，但预先明确地执行它通常更清晰/更透明）。 （我认为，因为您的数据集是平衡的，您可以首先在块结构上折叠/平均，然后用减少的数据集绘制整个图 - 它根本不应该改变线性回归的结果，至少不会平均值......无论如何，任何统计比较都应该在块级摘要上进行......）

df.sum <- (df.melted
    %>% group_by(Location,treatment,variable)
    %>% summarise(value=mean(value),yield_sd=sd(yield),
                  ## collapse yield to mean *after* computing sd!
                  yield=mean(yield))
)

剧情：

(ggplot(df.melted,
       aes(x = value, y = yield, color = treatment))
    + stat_smooth(method = "lm", col = "darkslategrey", se=FALSE)
    + stat_poly_eq(
          formula = y ~ x,
          ## aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
          aes(group=1, label =  ..rr.label..), 
          parse = TRUE)
    + theme_classic()
    + scale_shape(guide=FALSE)
    + geom_point(data=df.sum)
    + geom_errorbar(data=df.sum,
                    aes(ymax = yield+yield_sd, ymin = yield-yield_sd),
                    width = 0.05)
    + facet_wrap(~variable,scale="free_x")
)

（将group=1 添加到stat_poly_eq() 美学意味着我们每个方面只绘制一个R^2 值）

由于您不再将形状美学用于任何事物，您可以考虑使用它来显示 Location 变量...