0

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

我尝试使用此处列出的方法 从ggplot函数中提取信息，为 geom_density_ridges 添加平均值，并在geom_density_ridges 上画线，从中提取数据，但它没有提供我想要的所有信息 。stat_density_ridges()stat_density_ridgesggplot_buildpurrr::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 中的管道，因此如果解决方案中有这些，可能需要一些简单的解释。

4

1 回答 1

1

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

基本上，您可以group_by根据sample.number. 然后用于计算每个组内summarise列的汇总指标。x

要计算密度，您可以density()使用summarise. 这将返回一个包含(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)
于 2019-10-30T19:52:21.957 回答