R - Lattice xyplot - How do you add error bars to

Question

I'm posting this question because the very similar question here has not been answered until now.

I have been asked to plot the mean +/- SEM of my whole cohort of patients over the xyplot() that depicts the values of all patients. The data used represents intraoperative cardiovascular findings from patients undergoing surgery.

This is my data.frame called df

dput(df)
structure(list(Name = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 
5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L), .Label = c("DE", "JS", "KG", "MK", "TG", "WT"), class = "factor"), 
    Time = structure(c(1L, 2L, 3L, 4L, 5L, 7L, 8L, 1L, 2L, 3L, 
    4L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 
    6L, 7L, 8L, 2L, 3L, 4L, 5L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 
    7L, 8L), .Label = c("T1", "T2", "T3", "T4", "T5", "T6", "T7", 
    "T8"), class = "factor"), Dobut = structure(c(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, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L), .Label = c("No", "Yes"
    ), class = "factor"), DobutDose = c(NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    4L, 6L, 8L, 8L, 8L, 8L, 8L, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, 5L, 5L, NA), CI = c(1.4, 2.3, 1.3, 1.8, 2.1, 
    2, 2.1, 2.1, 2.3, 1.9, 1.6, 2, 2.4, 2.7, 2.6, 2.7, 2.6, 2.3, 
    2.4, 2.6, 0.9, 2.5, 2.1, 1.6, 1.5, 1.8, 2, 2, 1.9, 2.1, 2.3, 
    2, 2.4, 2.3, 2.6, 2.4, 2, 2.2, 1.6, 2.1, 2.5, 2.8), SvO2 = c(57L, 
    65L, 47L, 45L, 51L, 60L, 56L, 70L, 85L, 75L, 79L, 82L, 73L, 
    77L, 78L, 73L, 71L, 73L, 80L, 74L, 41L, 66L, 51L, 51L, 49L, 
    54L, 68L, 48L, 80L, 70L, 71L, 69L, 74L, 79L, 77L, 77L, 75L, 
    74L, 70L, 79L, 80L, 79L), SVRI = c(4000L, 1983L, 4000L, 2444L, 
    1981L, 2120L, 2514L, 2971L, 2157L, 3747L, 4300L, 3200L, 2867L, 
    1778L, 1169L, 1215L, 1262L, 1461L, 1600L, 1692L, 4978L, 1760L, 
    2019L, 2650L, 2827L, 2356L, 1800L, 2840L, 2063L, 2248L, 1948L, 
    2160L, 1733L, 2296L, 2677L, 2100L, 2640L, 2655L, 3950L, 2210L, 
    2848L, 2543L), MAP = c(80L, 65L, 86L, 74L, 67L, 65L, 74L, 
    90L, 70L, 90L, 96L, 94L, 100L, 82L, 60L, 61L, 62L, 62L, 69L, 
    71L, 70L, 71L, 77L, 73L, 75L, 77L, 61L, 85L, 65L, 74L, 70L, 
    67L, 69L, 74L, 92L, 71L, 88L, 93L, 89L, 79L, 97L, 97L), CVP = c(10L, 
    8L, 21L, 19L, 15L, 12L, 8L, 12L, 8L, 11L, 10L, 14L, 14L, 
    22L, 22L, 20L, 21L, 20L, 21L, 16L, 14L, 16L, 24L, 20L, 22L, 
    24L, 16L, 14L, 16L, 15L, 14L, 13L, 17L, 8L, 5L, 8L, 22L, 
    20L, 20L, 21L, 8L, 8L), PAP = c(23L, 22L, 36L, 36L, 34L, 
    32L, 22L, 33L, 28L, 36L, 36L, 40L, 37L, 37L, 40L, 35L, 35L, 
    34L, 38L, 36L, 45L, 43L, 55L, 49L, 52L, 54L, 43L, 47L, 27L, 
    25L, 23L, 22L, 28L, 21L, 20L, 25L, 33L, 33L, 38L, 35L, 33L, 
    29L), PCWP = c(15L, 11L, 28L, 26L, 23L, 21L, 11L, 26L, NA, 
    NA, 25L, 25L, NA, 27L, NA, NA, NA, NA, NA, NA, 30L, NA, NA, 
    NA, NA, NA, NA, NA, 19L, NA, NA, NA, NA, NA, 16L, NA, NA, 
    NA, NA, NA, NA, NA)), .Names = c("Name", "Time", "Dobut", 
"DobutDose", "CI", "SvO2", "SVRI", "MAP", "CVP", "PAP", "PCWP"
), class = "data.frame", row.names = c(NA, -42L))

Now the first xyplot I made for the variable CI looks like this

require(lattice)
xyplot(CI~Time, groups=Name, data=df, ty=c("l", "p"),
+        ,xlab="Measurement Time Point", 
ylab=expression("CI"~(l/min/m^"2")), main="Cardiac Index")

Now I was able to add the mean (black line) of the whole cohort, by doing the following

xyplot(CI~Time, groups=Name, data=df, ty=c("l", "p"), 
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.linejoin(x, y, horizontal = FALSE,..., col="black", lty=1, lwd=4)
       }
       ,xlab="Measurement Time Point", 
ylab=expression("CI"~(l/min/m^"2")), main="Cardiac Index")

Now I'd like to add +/- SE to the mean as a line above/below the mean, but nowhere can I find how to do this.

What I can do is using the latticeExtra package is add the loess line +/- SE, as below, but that's not the correct mathematical function I'm looking for. I've left the mean line in there to illustrate the difference between the two.

require(latticeExtra)
xyplot(CI~Time, groups=Name, data=df, ty=c("l", "p"), 
+        panel = function(x, y, ...) {
+            panel.xyplot(x, y, ...)
+            panel.linejoin(x, y, horizontal = FALSE,..., col="black", lty=1, lwd=4)
+            panel.smoother(x,y,se=TRUE, col.se="grey")
+        }
+        ,xlab="Measurement Time Point", 
ylab=expression("CI"~(l/min/m^"2")), main="Cardiac Index")

I have performed an extensive search through SO and the internet, but I haven't been able to find the right function to do this.

Help is very much appreciated! Thanks.

Answer 1

You could create your own panel function to plot a +/- SD region. For example

#new panel function
panel.se <- function(x, y, col.se=plot.line$col, alpha.se=.25, ...) {
    plot.line <- trellis.par.get("plot.line")
    xs <- if(is.factor(x)) {
       factor(c(levels(x) , rev(levels(x))), levels=levels(x))
    } else {
       xx <- sort(unique(x))
       c(xx, rev(xx))
    }
    means <- tapply(y,x, mean, na.rm=T)
    vars <- tapply(y,x, var, na.rm=T)
    Ns <- tapply(!is.na(y),x, sum)
    ses <- sqrt(vars/Ns)
    panel.polygon(xs, c(means+ses, rev(means-ses)), col=col.se, alpha=alpha.se)
}

and then you can use it like

#include new panel function
xyplot(CI~Time, groups=Name, data=df, ty=c("l", "p"), 
       panel = function(x, y, ...) {
           panel.se(x,y, col.se="grey")
           panel.xyplot(x, y, ...)
           panel.linejoin(x, y, horizontal = FALSE,..., col="black", lty=1, lwd=4)

       }
       ,xlab="Measurement Time Point", 
ylab=expression("CI"~(l/min/m^"2")), main="Cardiac Index")

which results in