2

我发布这个问题是因为直到现在还没有回答这里非常相似的问题。

我被要求将我的整个患者队列的平均值 +/- SEM 绘制在xyplot()描述所有患者值的图上。使用的数据代表接受手术的患者的术中心血管发现。

这是我的data.frame称呼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))

现在xyplot我为变量做的第一个CI看起来像这样

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

第一个 xyplot

现在我可以通过执行以下操作添加整个队列的平均值（黑线）

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

第二个xyplot

现在我想将 +/- SE 添加到平均值之上/之下的一条线，但我无法找到如何做到这一点。

我能做的是使用latticeExtra包添加黄土线 +/- SE，如下所示，但这不是我正在寻找的正确数学函数。我在那里留下了平均线来说明两者之间的区别。 

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

第三个xyplot

我已经通过 SO 和互联网进行了广泛的搜索，但我无法找到合适的功能来执行此操作。

非常感谢您的帮助！谢谢。

4

1 回答 1

3

您可以创建自己的面板函数来绘制 +/- SD 区域。例如

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

然后你可以像这样使用它

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

这导致

在此处输入图像描述

于 2014-10-13T05:47:09.977 回答