Drawing SE in xyplot with errorbars

Question

I am trying to construct a simple XY-Graph with the milk production (called FCM) of two different groups of cows (from the output I got from the mixed model, using the lsmeans and SE). I was able to construct the plot displaying the lsmeans using the xyplot function in lattice:

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

I now want to add the error bars. I tried some things with the panel.arrow function, just copying and pasting from other examples but didn´t get any further.

I would really appreciate some help!

My lsmeans2 dataset:

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

Answer 1

Here, a small panel function that allows plotting mid points, lower end points, and upper end points (based on acylam's solution). No "manual" scaling of the y-axis limits required.

Note the LHS of the formula in the xyplot() call.

pnf.mean.bars <- function(x, y, ...) {
  end <- length(y)
  end.mean <- end %/% 3L
  end.lower <- 2L * end.mean
  panel.xyplot(x, y[1:end.mean], ...)
  panel.arrows(x, y[1:end.mean], x, y[(end.mean + 1L):end.lower], length = 0.05, angle = 90)
  panel.arrows(x, y[1:end.mean], x, y[(end.lower + 1L):end], length = 0.05, angle = 90)
}

xyplot(c(lsmean, lower.CL, upper.CL) ~ Time.lab | Group, data = lsmeansSub,
       panel = pnf.mean.bars, type = "b", pch = 16, lwd = 2,
       scales = list(x = list(alternating = FALSE, rot = 90)),
       ylab = "FCM (kg/day)", xlab = "")

Answer 2

In case the two groups should be plotted in separate panels, this could be a solution.

It uses two of lattice's lesser known features, packet.number() and subscripts.

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

Answer 3

Is there a particular reason you want to use xplot? ggplot2 is much easier to work with and prettier. Here's an example of what I think you want.

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

Answer 4

For completeness, here is a solution using 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)))

The length argument in panel.arrows changes the width of the error heads. You can fiddle around with this parameter to get a width you like.

Notice that even though you had lsmeans2[order(lsmeans2$Time),] when specifying the data =, the ordering of Time is still wrong. This is because Time is a factor, and R doesn't know you want it to order by the numerical suffix of wk. This means, that it will sort wk10 before wk2, because 1 is smaller than 2. You can use this little trick below to order it correctly:

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

Note that even after reordering the the levels of "Time", I still need to use the sorted data for the data = argument. This is because xyplot plots the points in the order that appears in the dataset, not the order of the factor levels.