Рисование SE в xyplot с погрешностями

Я пытаюсь построить простой XY-график с производством молока (называемым FCM) двух разных групп коров (из результатов, которые я получил из смешанной модели, используя lsmeans и SE). Мне удалось построить график, отображающий lsmeans, используя функцию xyplot в решетке:

library(lattice)    
xyplot(lsmean~Time, type="b", group=Group, data=lsmeans2[order(lsmeans2$Time),],
       pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
       ylab="FCM (kg/day)", xlab="Week", lwd=2, 
       key=list(space="top",
                lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2),
                text=list(c("Confinement Group","Pasture Group"), cex=0.8)))

Теперь я хочу добавить полосы ошибок. Я пробовал кое-что с функцией panel.arrow, просто копировал и вставлял из других примеров, но дальше дело не пошло.

Я был бы очень признателен за помощь!

Мой набор данных lsmeans2:

Group Time lsmean SE df lower.CL upper.CL
Stall wk1  26.23299 0.6460481 59 24.19243 28.27356
Weide wk1  25.12652 0.6701080 58 23.00834 27.24471
Stall wk10 21.89950 0.6460589 59 19.85890 23.94010
Weide wk10 18.45845 0.6679617 58 16.34705 20.56986
Stall wk2  25.38004 0.6460168 59 23.33957 27.42050
Weide wk2  22.90409 0.6679617 58 20.79269 25.01549
Stall wk3  25.02474 0.6459262 59 22.98455 27.06492
Weide wk3  24.05886 0.6679436 58 21.94751 26.17020
Stall wk4  23.91630 0.6456643 59 21.87694 25.95565
Weide wk4  22.23608 0.6678912 58 20.12490 24.34726
Stall wk5  23.97382 0.6493483 59 21.92283 26.02481
Weide wk5  18.14550 0.6677398 58 16.03480 20.25620
Stall wk6  24.48899 0.6456643 59 22.44963 26.52834
Weide wk6  19.40022 0.6697394 58 17.28319 21.51724
Stall wk7  24.98107 0.6459262 59 22.94089 27.02126
Weide wk7  19.71200 0.6677398 58 17.60129 21.82270
Stall wk8  22.65167 0.6460168 59 20.61120 24.69214
Weide wk8  19.35759 0.6678912 58 17.24641 21.46877
Stall wk9  22.64381 0.6460481 59 20.60324 24.68438
Weide wk9  19.26869 0.6679436 58 17.15735 21.38004

r standard-error lattice

Melanie S 23.10.2016 источник

comment

Не могли бы вы показать dput(lsmeans2)? - Christoph 23.10.2016

Ответы (2)

arrow_upward
2
arrow_downward

Для полноты вот решение с использованием xyplot:

# Reproducible data
lsmeans2 = structure(list(Group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Stall", 
"Weide"), class = "factor"), Time = structure(c(1L, 1L, 2L, 2L, 
3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 
10L), .Label = c("wk1", "wk10", "wk2", "wk3", "wk4", "wk5", "wk6", 
"wk7", "wk8", "wk9"), class = "factor"), lsmean = c(26.23299, 
25.12652, 21.8995, 18.45845, 25.38004, 22.90409, 25.02474, 24.05886, 
23.9163, 22.23608, 23.97382, 18.1455, 24.48899, 19.40022, 24.98107, 
19.712, 22.65167, 19.35759, 22.64381, 19.26869), SE = c(0.6460481, 
0.670108, 0.6460589, 0.6679617, 0.6460168, 0.6679617, 0.6459262, 
0.6679436, 0.6456643, 0.6678912, 0.6493483, 0.6677398, 0.6456643, 
0.6697394, 0.6459262, 0.6677398, 0.6460168, 0.6678912, 0.6460481, 
0.6679436), df = c(59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 
58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L), lower.CL = c(24.19243, 
23.00834, 19.8589, 16.34705, 23.33957, 20.79269, 22.98455, 21.94751, 
21.87694, 20.1249, 21.92283, 16.0348, 22.44963, 17.28319, 22.94089, 
17.60129, 20.6112, 17.24641, 20.60324, 17.15735), upper.CL = c(28.27356, 
27.24471, 23.9401, 20.56986, 27.4205, 25.01549, 27.06492, 26.1702, 
25.95565, 24.34726, 26.02481, 20.2562, 26.52834, 21.51724, 27.02126, 
21.8227, 24.69214, 21.46877, 24.68438, 21.38004)), .Names = c("Group", 
"Time", "lsmean", "SE", "df", "lower.CL", "upper.CL"), class = "data.frame", row.names = c(NA, 
-20L))

xyplot(lsmean~Time, type="b", group=Group, data=lsmeans2[order(lsmeans2$Time),],
       panel = function(x, y, ...){
         panel.arrows(x, y, x, lsmeans2$upper.CL, length = 0.15,
                      angle = 90, col=c("darkorange","darkgreen"))
         panel.arrows(x, y, x, lsmeans2$lower.CL, length = 0.15,
                      angle = 90, col=c("darkorange","darkgreen"))
         panel.xyplot(x,y, ...)
       },
       pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
       ylab="FCM (kg/day)", xlab="Week", lwd=2, 
       key=list(space="top",
                lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2),
                text=list(c("Confinement Group","Pasture Group"), cex=0.8)))

Аргумент длины в panel.arrows изменяет ширину заголовков ошибок. Вы можете возиться с этим параметром, чтобы получить желаемую ширину.

Обратите внимание, что хотя у вас было lsmeans2[order(lsmeans2$Time),] при указании data =, порядок Time по-прежнему неправильный. Это связано с тем, что Time является фактором, и R не знает, что вы хотите, чтобы он упорядочивался по числовому суффиксу wk. Это означает, что он будет сортировать wk10 раньше, чем wk2, потому что 1 меньше 2. Вы можете использовать этот маленький трюк ниже, чтобы упорядочить его правильно:

# Order first by the character lenght, then by Time
Timelevels = levels(lsmeans2$Time) 
Timelevels = Timelevels[order(nchar(Timelevels), Timelevels)]

# Reorder the levels
lsmeans2$Time = factor(lsmeans2$Time, levels = Timelevels)

# Create Subset
lsmeansSub = lsmeans2[order(lsmeans2$Time),]

xyplot(lsmean~Time, type="b", group=Group, data=lsmeansSub,
       panel = function(x, y, yu, yl, ...){
         panel.arrows(x, y, x, lsmeansSub$upper.CL, length = 0.15,
                      angle = 90, col=c("darkorange","darkgreen"))
         panel.arrows(x, y, x, lsmeansSub$lower.CL, length = 0.15,
                      angle = 90, col=c("darkorange","darkgreen"))
         panel.xyplot(x, y, ...)
       },
       pch=16, ylim=c(10,35), col=c("darkorange","darkgreen"), 
       ylab="FCM (kg/day)", xlab="Week", lwd=2, 
       key=list(space="top",
                lines=list(col=c("darkorange","darkgreen"),lty=c(1,1),lwd=2),
                text=list(c("Confinement Group","Pasture Group"), cex=0.8)))

Обратите внимание, что даже после изменения порядка уровней «Время» мне все равно нужно использовать отсортированные данные для аргумента data =. Это связано с тем, что xyplot отображает точки в порядке, указанном в наборе данных, а не в порядке уровней факторов.

acylam 23.10.2016

comment

Большое спасибо! Я все еще учусь, и это мне очень помогло! - Melanie S; 25.10.2016

comment

@MelanieS Если вы считаете, что ответ полезен, примите его, чтобы другие могли его увидеть. - acylam; 31.08.2017

arrow_upward
1
arrow_downward

Есть ли особая причина, по которой вы хотите использовать xplot? ggplot2 с ним намного проще работать и он красивее. Вот пример того, что я думаю, вы хотите.

#load ggplot2
library(ggplot2)

#load data
d = structure(list(Group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Stall", 
"Weide"), class = "factor"), Time = structure(c(1L, 1L, 2L, 2L, 
3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 
10L), .Label = c("wk1", "wk10", "wk2", "wk3", "wk4", "wk5", "wk6", 
"wk7", "wk8", "wk9"), class = "factor"), lsmean = c(26.23299, 
25.12652, 21.8995, 18.45845, 25.38004, 22.90409, 25.02474, 24.05886, 
23.9163, 22.23608, 23.97382, 18.1455, 24.48899, 19.40022, 24.98107, 
19.712, 22.65167, 19.35759, 22.64381, 19.26869), SE = c(0.6460481, 
0.670108, 0.6460589, 0.6679617, 0.6460168, 0.6679617, 0.6459262, 
0.6679436, 0.6456643, 0.6678912, 0.6493483, 0.6677398, 0.6456643, 
0.6697394, 0.6459262, 0.6677398, 0.6460168, 0.6678912, 0.6460481, 
0.6679436), df = c(59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 
58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L, 59L, 58L), lower.CL = c(24.19243, 
23.00834, 19.8589, 16.34705, 23.33957, 20.79269, 22.98455, 21.94751, 
21.87694, 20.1249, 21.92283, 16.0348, 22.44963, 17.28319, 22.94089, 
17.60129, 20.6112, 17.24641, 20.60324, 17.15735), upper.CL = c(28.27356, 
27.24471, 23.9401, 20.56986, 27.4205, 25.01549, 27.06492, 26.1702, 
25.95565, 24.34726, 26.02481, 20.2562, 26.52834, 21.51724, 27.02126, 
21.8227, 24.69214, 21.46877, 24.68438, 21.38004)), .Names = c("Group", 
"Time", "lsmean", "SE", "df", "lower.CL", "upper.CL"), class = "data.frame", row.names = c(NA, 
-20L))

#fix week
library(stringr)
library(magrittr)
d$Time %<>% as.character() %>% str_replace(pattern = "wk", replacement = "") %>% as.numeric()

#plot
ggplot(d, aes(Time, lsmean, color = Group, group = Group)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2) +
  geom_line() +
  ylim(10, 35) +
  scale_x_continuous(name = "Week", breaks = 1:10) +
  ylab("FCM (kg/day)") +
  scale_color_discrete(label = c("Confinement Group","Pasture Group"))

CoderGuy123 23.10.2016

comment

Большое спасибо! Однако время не в правильном порядке, но userR решил это в своем примере. Я очень благодарна за быстрый ответ! - Melanie S; 25.10.2016

comment

Я обновил код, чтобы исправить порядок недели. Обратите внимание, что я только что использовал данные, которые вы мне дали, и ваши данные были не в порядке. Я удалил wk из строк, потому что это было лишним; единица времени уже обозначена меткой X. - CoderGuy123; 25.10.2016

Рисование SE в xyplot с погрешностями

Ответы (2)

Вопросы по теме