График линейного регрессионного анализа с погрешностью для изменчивости - PullRequest
Новая
       21

График линейного регрессионного анализа с погрешностью для изменчивости

1 голос
/ 30 октября 2019

Я хотел сделать графики, похожие на рисунок 1 (источник: ссылка )

enter image description here

На рисунке 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")

1 Ответ

2 голосов
/ 30 октября 2019 
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 на фасет)

enter image description here

Поскольку вы больше не используете эстетику формы для чего-либо, вы можете использовать его, чтобы показать переменную Location ...

...