r - 如果可能，在 R - 基 R 中着色置信区间

Question

我正在比较使用 LOESS 回归的两条线。我想清楚地显示两条线的置信区间，我遇到了一些困难。

我尝试过使用各种线型和颜色，但在我看来结果仍然是忙碌和混乱。我认为置信区间之间的阴影可能会使事情变得更清晰，但是考虑到我的编码到目前为止的结构，我在解决这个问题时遇到了一些困难。我已经包含了生成的图、Analysis5k 和 Analysis5kz 两组的数据，以及到目前为止的代码。

我已经看到了一些示例，其中两个多边形重叠以显示置信区间重叠的位置，这似乎是呈现数据的好方法。如果有一种方法可以在两个置信区间共享的区域中绘制多边形，那可能是呈现数据的另一种好方法。

我了解应该如何完成多边形的基本概念，但我发现的示例已应用于更简单的线条和数据。到目前为止，对于一些糟糕的组织来说，部分原因是我自己的错，但由于这一步基本上是我数据呈现的最后润色，我真的不想从头开始重新工作。

非常感谢任何帮助或见解。

更新

我更新了标题。我收到了一些使用 ggplot 的好例子，虽然我想在未来使用 ggplot，但到目前为止我只处理了基础 R。对于这个特定的项目，如果可能的话，想尝试将其保留在基础 R 中。

分析5k

Period  15p5    Total_5plus
-4350   0.100529101 12.6
-3900   0.4 20
-3650   0.0625  9.6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3650   0.174757282 20.6
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-3400   0.208333333 24
-4233   0.184027778 19.2
-3650   0.285714286 12.6
-4350   0.166666667 6

分析5kz

Period  15p5    Total_5plus
-4350   0.100529101 12.6
-4350   0   5
-3900   0.4 20
-3650   0.0625  9.6
-3400   0   6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3650   0   28
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0   25
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3500   0   5
-3650   0.174757282 20.6
-3800   0   9
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-4117   0   8
-4350   0   8
-3400   0.208333333 24
-4233   0.184027778 19.2
-3025   0   7
-3650   0.285714286 12.6
-4350   0.166666667 6

代码

  ppi <- 300 
  png("5+ KC shaded CI.png", width=6*ppi, height=6*ppi, res=ppi) 
  library(Hmisc) 
  Analysis5k <- read.csv(file.choose(), header = T) 
  Analysis5kz <- read.csv(file.choose(), header = T)
  par(mfrow = c(1,1), pty = "s", oma=c(1,2,1,1), mar=c(4,4,2,2)) 
  plot(X15p5 ~ Period, Analysis5kz, xaxt = "n", yaxt= "n", ylim=c(-0.2,0.7), xlim=c(-5000,-2500), xlab = "Years B.P.", ylab = expression(''[15]*'p'[5]), main = "") 
  vx <- seq(-5000,-2000, by = 500) 
  vy <- seq(-0.2,0.7, by = 0.1) 
  axis(1, at = vx) 
  axis(2, at = vy) 
  a5k <- order(Analysis5k$Period) 
  a5kz <- order(Analysis5kz$Period)
  Analysis5k.lo <- loess(X15p5 ~ Period, Analysis5k, weights = Total_5plus, span = 0.6) 
  Analysis5kz.lo <- loess(X15p5 ~ Period, Analysis5kz, weights = Total_5plus, span = 0.6)      
  pred5k <- predict(Analysis5k.lo, se = TRUE) 
  pred5kz <- predict(Analysis5kz.lo, se = TRUE)      
  lines(Analysis5k$Period[a5k], pred5k$fit[a5k], col="blue", lwd=2) 
  lines(Analysis5kz$Period[a5kz], pred5kz$fit[a5kz], col="skyblue", lwd=2)          
  lines(Analysis5K$Period[a5K], pred5K$fit[a5K] - qt(0.975, pred5K$df)*pred5K$se[a5K],col="blue",lty=2) 
  lines(Analysis5K$Period[a5K], pred5K$fit[a5K] + qt(0.975, pred5K$df)*pred5K$se[a5K],col="blue",lty=2)      
  lines(Analysis5Kz$Period[a5Kz], pred5Kz$fit[a5Kz] - qt(0.975, pred5Kz$df)*pred5Kz$se[a5Kz],col="skyblue",lty=2) 
  lines(Analysis5Kz$Period[a5Kz], pred5Kz$fit[a5Kz] + qt(0.975, pred5Kz$df)*pred5Kz$se[a5Kz],col="skyblue",lty=2)
  abline(h=0.173, lty=3) 
  abline(v=-4700, lty=3)
  abline(v=-4000, lty=3)
  abline(v=-3000, lty=3)
  minor.tick(nx=5, ny=4, tick.ratio=0.5) 
  dev.off()

score 5 · Accepted Answer

这是基于您的代码的基本图的解决方案。

诀窍polygon是您必须在一个向量中提供 2 倍的 x 坐标，一次按正常顺序，一次按相反顺序（使用 function rev），并且您必须提供 y 坐标作为上限的向量，然后是下限相反的顺序。

我们使用该adjustcolor功能使标准颜色透明。

library(Hmisc) 
ppi <- 300 
par(mfrow = c(1,1), pty = "s", oma=c(1,2,1,1), mar=c(4,4,2,2)) 
plot(X15p5 ~ Period, Analysis5kz, xaxt = "n", yaxt= "n", ylim=c(-0.2,0.7), xlim=c(-5000,-2500), xlab = "Years B.P.", ylab = expression(''[15]*'p'[5]), main = "") 
vx <- seq(-5000,-2000, by = 500) 
vy <- seq(-0.2,0.7, by = 0.1) 
axis(1, at = vx) 
axis(2, at = vy) 
a5k <- order(Analysis5k$Period) 
a5kz <- order(Analysis5kz$Period)
Analysis5k.lo <- loess(X15p5 ~ Period, Analysis5k, weights = Total_5plus, span = 0.6) 
Analysis5kz.lo <- loess(X15p5 ~ Period, Analysis5kz, weights = Total_5plus, span = 0.6)      
pred5k <- predict(Analysis5k.lo, se = TRUE) 
pred5kz <- predict(Analysis5kz.lo, se = TRUE)      

polygon(x = c(Analysis5k$Period[a5k], rev(Analysis5k$Period[a5k])),
        y = c(pred5k$fit[a5k] - qt(0.975, pred5k$df)*pred5k$se[a5k], 
              rev(pred5k$fit[a5k] + qt(0.975, pred5k$df)*pred5k$se[a5k])),
        col =  adjustcolor("dodgerblue", alpha.f = 0.10), border = NA)

polygon(x = c(Analysis5kz$Period[a5kz], rev(Analysis5kz$Period[a5kz])),
        y = c(pred5kz$fit[a5kz] - qt(0.975, pred5kz$df)*pred5kz$se[a5kz], 
              rev( pred5kz$fit[a5kz] + qt(0.975, pred5kz$df)*pred5kz$se[a5kz])),
        col =  adjustcolor("orangered", alpha.f = 0.10), border = NA)

lines(Analysis5k$Period[a5k], pred5k$fit[a5k], col="dodgerblue", lwd=2) 
lines(Analysis5kz$Period[a5kz], pred5kz$fit[a5kz], col="orangered", lwd=2)   

abline(h=0.173, lty=3) 
abline(v=-4700, lty=3)
abline(v=-4000, lty=3)
abline(v=-3000, lty=3)
minor.tick(nx=5, ny=4, tick.ratio=0.5)

score 2 · Accepted Answer

这是使用 ggplot 的一种方法：

(1) 对两个数据集应用黄土平滑

library(dplyr)
df.lo <- lapply(datlist, function(x)loess(X15p5 ~ Period, data=x, weights = Total_5plus, span = 0.6))

(2) 创建一个新的 data.frame，扩展 data.set 的 min (-4350) 和 max Period (-2715)：

nd1 <- nd2 <- expand.grid(Period=seq(-4350, -2715, length=100))

(3) 预测每个 loess 平滑器的 fit 和 se 并绑定到单个 data.frame 中：

nd1[,c("fit", "se")] <- predict(df1.lo[[1]], newdata=nd1, se=T)[1:2]
nd1 <- nd1 %>% mutate(group="5k")
nd2[,c("fit", "se")] <- predict(df2.lo[[2]], newdata=nd1, se=T)[1:2]
nd2 <- nd2 %>% mutate(group="5kz")

ndata <- rbind(nd1, nd2)

(4) 用预测的数据，ggplot2::geom_ribbon用来显示重叠的se：

library(ggplot2)
p <- ggplot(ndata, aes(Period, fit)) + 
  geom_line(aes(colour=group)) + 
  geom_ribbon(aes(ymin=fit-1.96*se, ymax=fit+1.96*se, fill=group), alpha=.2) 

p

(5)添加数据点和abline：

dat <- do.call(rbind, datlist)
p + 
  geom_point(data=dat, aes(y=X15p5, shape=as.factor(group)), alpha=.2) + 
  geom_hline(yintercept=0.173, linetype="dotted") + 
  geom_vline(xintercept=c(-4700, -4000, -3000), linetype="dotted") +
  ylab("X15p5") + 
  theme_bw()

源数据datlist是两个data.frames“Analysis5k”和“Analysis5kz”的列表。输出如下：

structure(list(`5k` = structure(list(Period = c(-4350L, -3900L, 
-3650L, -3900L, -3958L, -4350L, -3400L, -3650L, -3400L, -3500L, 
-3600L, -3808L, -3900L, -3958L, -4058L, -3500L, -4117L, -4350L, 
-4350L, -3650L, -2715L, -4350L, -3500L, -4350L, -3650L, -2715L, 
-3500L, -3650L, -3400L, -4233L, -3650L, -4350L), X15p5 = c(0.100529101, 
0.4, 0.0625, 0.126984127, 0.133333333, 0.150943396, 0.146341463, 
0.255319149, 0.222222222, 0.245014245, 0.125, 0.1, 0.160493827, 
0.238095238, 0.2, 0.086956522, 0.141414141, 0.171038825, 0.166666667, 
0.143798024, 0.137931034, 0.235588972, 0.228840125, 0.041666667, 
0.174757282, 0.377777778, 0.2, 0.078947368, 0.208333333, 0.184027778, 
0.285714286, 0.166666667), Total_5plus = c(12.6, 20, 9.6, 16.8, 
5, 10.6, 8.2, 9.4, 9, 39, 8, 20, 18, 7, 5, 28.75, 6.6, 31.76666667, 
6, 30.36666667, 7.25, 26.6, 79.75, 8, 20.6, 11.25, 7.5, 7.6, 
24, 19.2, 12.6, 6), group = c("5k", "5k", "5k", "5k", "5k", "5k", 
"5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", 
"5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", "5k", 
"5k", "5k", "5k", "5k")), .Names = c("Period", "X15p5", "Total_5plus", 
"group"), row.names = c(NA, 32L), class = "data.frame"), `5kz` = 
structure(list(
    Period = c(-4350L, -4350L, -3900L, -3650L, -3400L, -3900L, 
    -3958L, -4350L, -3400L, -3650L, -3400L, -3500L, -3600L, -3650L, 
    -3808L, -3900L, -3958L, -4058L, -3500L, -3500L, -4117L, -4350L, 
    -4350L, -3650L, -2715L, -4350L, -3500L, -4350L, -3500L, -3650L, 
    -3800L, -2715L, -3500L, -3650L, -4117L, -4350L, -3400L, -4233L, 
    -3025L, -3650L, -4350L), X15p5 = c(0.100529101, 0, 0.4, 0.0625, 
    0, 0.126984127, 0.133333333, 0.150943396, 0.146341463, 0.255319149, 
    0.222222222, 0.245014245, 0.125, 0, 0.1, 0.160493827, 0.238095238, 
    0.2, 0, 0.086956522, 0.141414141, 0.171038825, 0.166666667, 
    0.143798024, 0.137931034, 0.235588972, 0.228840125, 0.041666667, 
    0, 0.174757282, 0, 0.377777778, 0.2, 0.078947368, 0, 0, 0.208333333, 
    0.184027778, 0, 0.285714286, 0.166666667), Total_5plus = c(12.6, 
    5, 20, 9.6, 6, 16.8, 5, 10.6, 8.2, 9.4, 9, 39, 8, 28, 20, 
    18, 7, 5, 25, 28.75, 6.6, 31.76666667, 6, 30.36666667, 7.25, 
    26.6, 79.75, 8, 5, 20.6, 9, 11.25, 7.5, 7.6, 8, 8, 24, 19.2, 
    7, 12.6, 6), group = c("5kz", "5kz", "5kz", "5kz", "5kz", 
    "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", 
    "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", 
    "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", 
    "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz", "5kz"
    )), .Names = c("Period", "X15p5", "Total_5plus", "group"), row.names = 33:73, class = "data.frame")), .Names = c("5k", 
"5kz"))

score 2 · Accepted Answer

我会提出一个 tidyverse 解决方案。在这种方法中，您首先创建一个函数来计算和提取所需的统计数据。nest然后，您使用、map该列表上的函数和结果创建一个列表列unnest。

您可以在http://r4ds.had.co.nz/many-models.html阅读有关此方法的更多信息。

library(tidyverse)

# create function to retrieve fit and se
pred_fun <- function(df) {
  model <- loess(`15p5` ~ Period, df, weights = Total_5plus, span = .6)
  preds <- predict(model, se = T)

  data_frame(fit = preds[["fit"]],
             se = preds[["se.fit"]])
}

# nest, map and unnest fits
nested <- bind_rows(df_5k, df_5kz) %>% 
  group_by(origin) %>% 
  nest() %>% 
  mutate(preds = map(data, pred_fun)) %>% 
  unnest(data, preds)


# plot result
ggplot(nested, aes(Period, `15p5`)) +
  geom_ribbon(aes(ymin = fit - 1.96 * se, ymax = fit + 1.96 * se, fill = origin),
              alpha = .2) +
  geom_point() +
  geom_line(aes(y = fit, colour = origin)) +
  scale_y_continuous(expand = c(.3, 0)) +
  scale_x_continuous(expand = c(.3, 0), breaks = scales::pretty_breaks(6)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  labs(x = "Years B.P.", y = expression(''[15]*'p'[5]))

当然，您可以编辑组的颜色，例如：

cols <- c(df_5k = "blue", df_5kz = "skyblue")

ggplot...
...
scale_fill_manual(values = cols) +
scale_color_manual(values = cols)

编辑：

由于我不知道如何使用基本图形执行您想要的操作，因此我会尝试使绘图看起来像基础，使用ggthemes::theme_base和更改点类型如下：

ggplot(nested, aes(Period, `15p5`)) +
  ggthemes::theme_base() +
  geom_hline(yintercept = 0.173, linetype = "dotted") +
  geom_vline(xintercept = c(-4700, -4000, -3000), linetype = "dotted") +
  geom_ribbon(aes(ymin = fit - 1.96 * se, ymax = fit + 1.96 * se, fill = origin),
              alpha = .2) +
  geom_point(shape = 1) +
  geom_line(aes(y = fit, colour = origin)) +
  scale_y_continuous(expand = c(.3, 0)) +
  scale_x_continuous(expand = c(.3, 0), breaks = scales::pretty_breaks(6)) +
  theme(legend.position = "bottom") +
  labs(x = "Years B.P.", y = expression(''[15]*'p'[5]),
       colour = NULL, fill = NULL)

数据导入

df_5k <- "Period  15p5    Total_5plus
-4350   0.100529101 12.6
-3900   0.4 20
-3650   0.0625  9.6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3650   0.174757282 20.6
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-3400   0.208333333 24
-4233   0.184027778 19.2
-3650   0.285714286 12.6
-4350   0.166666667 6"

df_5k <- read_table2(df_5k) %>% 
  mutate(origin = "df_5k")

df_5kz <- "Period  15p5    Total_5plus
-4350   0.100529101 12.6
-4350   0   5
-3900   0.4 20
-3650   0.0625  9.6
-3400   0   6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3650   0   28
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0   25
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3500   0   5
-3650   0.174757282 20.6
-3800   0   9
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-4117   0   8
-4350   0   8
-3400   0.208333333 24
-4233   0.184027778 19.2
-3025   0   7
-3650   0.285714286 12.6
-4350   0.166666667 6"

df_5kz <- read_table2(df_5kz) %>% 
  mutate(origin = "df_5kz")

r - 如果可能，在 R - 基 R 中着色置信区间

3 回答 3

编辑：

数据导入

Related

Reference