geom_boxplot 样本 ID 中的异常值形状

geom_boxplot outlier shape from Sample ID

如何修改 geom_boxplot 中离群值的形状以随着时间的推移匹配样本 ID。 想象一下我有这种数据(这只是虚拟数据,代码可能不漂亮,但这就是我想出的):

# create dummy data
df <- data.frame()
set.seed(42)
os <- 0
sam <- 1
for (time in as.factor(c('T0', 'T1'))) {
  if (time == 'T1') {
    sam <- 1
  }
  for (group in as.factor(c('A','B'))) {
    for (pat in 1:10) {
      df[pat + os, 'Sample'] <- paste('P', pat, '_', sam, sep = '')
      df[pat + os, 'Time'] <- time
      df[pat + os, 'Group'] <- group
      df[pat + os, 'Value'] <- rnorm(1) + os
      # add outlier, they are the same in each group in this example,
      # but can differ in the real data set
      if (pat == 2 | pat == 9) {
        print(pat)
        df[pat + os, 'Value'] <- df[pat + os, 'Value'] + 10
      }
      sam <- sam + 1
    }
    os <- os + 10
  }
}

# mark outliers in table
df = df %>% 
  group_by(Group,Time)  %>%
  mutate(is_outlier = case_when(Value > quantile(Value)[4] + 1.5*IQR(Value) ~ TRUE,
                                Value < quantile(Value)[2] - 1.5*IQR(Value) ~ TRUE,
                                TRUE ~ FALSE))

这导致以下情节:

ggplot(df, aes(x = Time,
               y = Value,
               label = Time)) +
  geom_boxplot(outlier.colour =  'red',
               outlier.shape = 1,
               outlier.size = 2
  ) +
  facet_grid(~factor(Group),
             switch = 'x',
             scales = 'free_y')

目标:

我想要的是对于每个组 AB 我可以查看异常值是否相同。因此,例如 A T0 中显示的离群值与 A T1 中的相同。更具体地说,在 A T0 中被视为圆圈的异常值在 A T1 中应该是一个圆圈,而在 A T1 中的第二个异常值应该是任何其他形状(例如三角形)。由于我的原始数据有大约 5/6 个时间点,因此很高兴通过查看绘图来了解异常值是否仍然是异常值。 在某些情况下,我的原始数据集有大约 5-8 个异常值。

在组 B 中,我们可以重复使用与组 A 中相同的形状,尽管我们的样本 ID 与组 A 中不同。

我想使用三角形、圆形、Asterix 等基本形状(我知道形状有限,但对于我的数据集来说应该足够了)。我也知道我可以标记数据点,但我不想这样做。 不同的颜色也可以,但我更喜欢不同的形状。

我想我必须单独计算离群值,然后可能使用 geom_point 和 aes(shape = df$Sample) 之类的东西。但是我想不通。

是否有人根据我的虚拟数据提供了提示或解决方案? 那太棒了:-)

最佳TMC

我想出了一个非常丑陋的解决方案。我很确定有一种更漂亮的方法可以做到这一点,但这是完整的代码:

首先我们创建虚拟数据:

# start with an clean environment
rm(list=ls())  
# create a function to load or install all necessary libraries 
install.load.package <- function(x) {
  if (!require(x, character.only = TRUE))
    install.packages(x)
  require(x, character.only = TRUE)
}
package_vec <- c("ggplot2",
                 "dplyr"
)
sapply(package_vec, install.load.package)  

# now to the data
df <- data.frame()
set.seed(42)
os <- 0
sam <- 1
for (time in as.factor(c('T0', 'T1'))) {
  if (time == 'T1') {
    sam <- 1
  }
  for (group in as.factor(c('A','B'))) {
    for (pat in 1:10) {
      df[pat + os, 'Sample'] <- paste('P', pat, '_', sam, sep = '')
      df[pat + os, 'Time'] <- time
      df[pat + os, 'Group'] <- group
      df[pat + os, 'Value'] <- rnorm(1) + os
      # add outlier, they are the same in each group in this example,
      # but can differ in the real data set
      if (pat == 2 | pat == 9) {
        print(pat)
        df[pat + os, 'Value'] <- df[pat + os, 'Value'] + 10
      }
      sam <- sam + 1
    }
    os <- os + 10
  }
}

然后我们计算异常值如下,并创建一个新列,其中放置异常值的ID。如果不是异常值,则插入 'X'

# calculate outliers
df = df %>% 
  group_by(Group,Time)  %>%
  mutate(is_outlier = case_when(Value > quantile(Value)[4] + 1.5*IQR(Value) ~ as.character(Sample),
                                Value < quantile(Value)[2] - 1.5*IQR(Value) ~ as.character(Sample),
                                TRUE ~ as.character('X')))
df$Group <- as.factor(df$Group)

现在,我们将样本 ID 替换为数字。第一个异常值对获得数字 1,第二个异常值对获得数字 2,依此类推。如果离群值多于可用的“geom_points”形状,则必须调整代码。但是让我们假设我们没有超过 23 个异常值(我认为这是最大数量)。

for (group in levels(df$Group)) {
  count <- 1
  for (id in levels(as.factor(df$is_outlier[which(df$Group == group)]))) {
    if (id == 'X') {
      df[which(df$is_outlier == id), 'is_outlier'] <- as.character(NA)
    } else {
      df[which(df$is_outlier == id), 'is_outlier'] <- as.character(count)
      count <- count + 1
    }
  }
}

这会覆盖之前创建的列。它为 X 值引入 NA

现在我们可以绘制数据了

  ggplot(df, aes(x = Time,
              y = Value,
              label = Time)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(data = df,
             shape= as.numeric(df$is_outlier),
             color = 'red') +
  facet_grid(~factor(Group),
             switch = 'x',
             scales = 'free_y')

这导致了这个情节:

现在我们可以查看异常值是否在 T0T1 之间保持异常值。请注意,在组 B 中,我们使用相同的形状。但这些是完全不同的样本。必须修改绘图代码上方的代码来解决这一问题。但这样一来,我们可用的形状可能会更少。

如果你们有更流畅、更优雅的解决方案,我很乐意学习。

最佳TMC