如何一次从 R 中的核密度函数中提取多个样本的数据

Question

我有一个非常大的数据文件（>300k 行），每一行都是唯一样本（>3000 个样本）的一部分。我想为每个单独的样本生成一个核密度估计器，并将相关信息（最小值、最大值、密度估计器的最大概率、密度估计器的中值、密度估计器的平均值）与样本名称一起提取到一个单独的表中。

我尝试使用此处列出的方法从ggplot 函数stat_density_ridges() 中提取信息 Adding a mean to geom_density_ridges 和这里的draw line on geom_density_ridges 从stat_density_ridges 和ggplot_build 和purrr::pluck 中提取数据，但它没有提供我想要的所有信息。

以下生成一些类似于我想要的合成数据：

set.seed(1)
x = runif( 50, max = 40, min = 20 )
set.seed(2)
y = runif( 50, max = 300, min = 100 )
sample.number = c( rep( 1, 20 ), rep( 2, 15 ), rep( 3, 5 ), rep( 4, 10 ) )
d <- data.frame( x, y , sample.number ) 

ggplot 中显示分布的图：

ggplot( data = d, aes( x = x, y = as.factor( samples ) ) ) +
  labs( x = expression( paste( "x" ) ), 
    y = expression( paste( "sample number" ) ) ) +
  stat_density_ridges() 

我想最终得到一个包含以下信息的数据表： sample.name，max(x)，min(x)，内核密度估计器的最大高度是x位置，内核密度估计器的中值高度是x位置等。

我唯一能想到的就是创建一个漫长而艰巨的循环

sample.numbers <- rep( NA, times = max( d$sample.number ) )
max.x <- rep( NA, times = max( d$sample.number ) )
min.x <- rep( NA, times = max( d$sample.number ) )

for( i in 1:max( d$sample.number ) ) {
  temp.d = d[ d$sample.number == i, ]
  sample.numbers[ i ] = i
  max.x[ i ] = max( temp.d$x )
  min.x[ i ] = min( temp.d$x )
}

然后以某种方式添加一些创建密度估计器并从中提取信息的位。我猜 R 中的索引为我在使用 group_by 时拥有的数千个样本提供了一种更简单的方法来解决这个问题，但我无法弄清楚。请注意，我仍然无法理解 R 中的管道，因此如果解决方案中有这些，可能需要一些简单的解释。

Answer 1

有不同的方法可以做到这一点。在我看来，使用 dplyr 和管道运算符是最简单的方法。我尝试在代码中添加 cmets 以使其更易于理解。看看this dplyr cheat sheet。

基本上，您使用group_by 根据sample.number 将数据框分组。然后使用summarise 计算每个组内x 列的汇总指标。

要计算密度，您可以使用summarise 内的基数 R 中的density()。这将返回一个列表，其中包含密度函数的(x,y) 值样本。要从此密度函数中提取分位数，可以使用包spatstat。

观察：density() 计算取决于数据集的带宽值。由于我们将不同的组分开，因此每个组最终可能具有不同的带宽值。我使用函数bw.nrd 使用完整数据集估计带宽的单个值。然后我将这个单一带宽值用于所有计算。

# needed to extract quantile from a pdf computed with density()
library(spatstat)
# packages for data wrangling
library(plyr)
library(dplyr)
# ploting
library(ggplot2)
library(ggridges)

# creata data set
set.seed(1)
x = runif( 50, max = 40, min = 20 )
set.seed(2)
y = runif( 50, max = 300, min = 100 )
sample.number = c( rep( 1, 20 ), rep( 2, 15 ), rep( 3, 5 ), rep( 4, 10 ) )
d <- data.frame( x, y , sample.number )

# first compute bandwidth over all samples
# if you don't do this, each pdf in the table will have a different bandwidth
# bw.nrd is a function that computes bandwidth for a kernel density using a "rule of thumb" formula
# there are other functions that you can use to estimate bw
bw <- bw.nrd(d$x)

# create the table using the pipe operator and dplyr
# the pipe operator '%>%' takes what is on the left side and puts inside the function
# on the right side as an argument
d %>% 
  # group rows of 'd' by sample number (this is equivalent to your for loop)
  group_by(sample.number) %>%
  # before computing the summaries for each group, create a new column with the 
  # number of elements in each sample (the resulting DF still has 50 rows)
  mutate(n=n()) %>%
  # now remove rows that belong to groups with less than 5 elements (you can change the threshold value here)
  filter(n > 5) %>%
  # for each group in 'd' compute these summary metrics
  summarise(max.x=max(x),
            min.x=min(x), 
            max.density=max(density(x, bw = bw)$y),
            x.mode=density(x, bw = bw)$x[which(density(x, bw = bw)$y == max.density)],
            x.median=quantile(density(x, bw = bw), 0.5),
            median.density=density(x, bw = bw)$y[which(density(x, bw = bw)$x == x.median)])

# OUTPUT (note that sample.number == 3 was removed from the table)
#># A tibble: 3 x 7
#>  sample.number max.x min.x max.density x.mode x.median median.density
#>          <dbl> <dbl> <dbl>       <dbl>  <dbl>    <dbl>          <dbl>
#>1             1  39.8  21.2      0.0568   34.3     31.4         0.0503
#>2             2  38.7  20.3      0.0653   26.9     28.4         0.0628
#>3             4  36.4  20.5      0.0965   33.9     33.0         0.0939
#

# see the pdfs using stat_density_ridges
# (note that i am fixing the bandwidth)
ggplot( data = d, aes( x = x, y = as.factor( sample.number ) ) ) +
  labs( x = expression( paste( "x" ) ), 
        y = expression( paste( "sample number" ) ) ) +
  stat_density_ridges(bandwidth = bw)