在 xyplot 中用误差线绘制 SE

问题描述 投票：0回答：3

我正在尝试构建一个简单的 XY 图，其中包含两组不同奶牛的产奶量（称为 FCM）（根据我从混合模型获得的输出，使用 lsmeans 和 SE）。 我能够使用点阵中的 xyplot 函数构建显示 lsmeans 的图：

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 lattice standard-error
3个回答
2
投票

为了完整起见，这里是一个使用 

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 =
，时间的顺序仍然是错误的。这是因为时间是一个因素，而 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
按照数据集中出现的顺序绘制点，而不是因子水平的顺序。

1
投票

您想使用 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"))

0
投票

如果这两个组应该绘制在单独的面板中，这可能是一个解决方案。

它使用了lattice的两个鲜为人知的功能，packet.number()和下标。

##  acylam's part.
##  ==============

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

##  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),]






##  My part.
##  ========

##  "Beautify" labels of 'Time' (for plotting only; include spaces).
##  ----------------------------------------------------------------
lsmeansSub$Time.lab <- lsmeansSub$Time
levels(lsmeansSub$Time.lab) <- c("Week 1", "Week 2", "Week 3", "Week 4", "Week 5", "Week 6", "Week 7", "Week 8", "Week 9", "Week 10")

##  Panel function, using "packet.number()" for the two panels, and
##  "subscripts" for selecting the appropriate CI limits.
##  Note the additional arguments "limits" and "colors" (whilst
##  "subscripts" is from lattice).
##  -----------------------------------------------------------------
pnf <- function(x, y, limits, colors, subscripts, ...) {
    panel.arrows(x, y, x, limits[subscripts, "upper.CL"], length = 0.15, angle = 90, col = colors[packet.number()])
    panel.arrows(x, y, x, limits[subscripts, "lower.CL"], length = 0.15, angle = 90, col = colors[packet.number()])
    panel.xyplot(x, y, col = colors[packet.number()], ...)
}


##  Plot.
##  Note the passing of the specific values to "limits" and "colors".
##  -----------------------------------------------------------------
xyplot(lsmean ~ Time.lab | Group, data = lsmeansSub,
       limits = lsmeansSub[, c("lower.CL", "upper.CL")],
       colors = c("darkorange", "darkgreen"),
       panel = pnf, type = "b", pch = 16, lwd = 2,
       scales = list(x = list(alternating = FALSE, rot = 90)),
       ylim = c(10, 35), ylab = "FCM (kg/day)", xlab = "",
       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)))
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.