如何一次从 R 中的核密度函数中提取多个样本的数据
How can I extract data from a kernel density function in R for many samples at one time
我有一个非常大的数据文件（>300k 行），每一行都是一个独特样本（>3000 个样本）的一部分。我想为每个单独的样本生成一个核密度估计器，并将相关信息（最小值、最大值、密度估计器的最大概率、密度估计器的中值、密度估计器的平均值）提取到单独的 table 以及样本名称。
我已尝试使用此处列出的方法从 ggplot 函数 stat_density_ridges() 中提取信息
Adding a mean to geom_density_ridges and here 使用 purrr::pluck 从 stat_density_ridges 和 ggplot_build 中提取数据，但它没有提供我想要的所有信息。
下面生成一些类似于我想要的合成数据：
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()
我想得到一个包含以下信息的数据 table：
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 中的管道，因此如果解决方案中包含管道，则可能需要一些简单的解释。
有多种方法可以做到这一点。在我看来，使用 dplyr 和管道运算符是最简单的方法。我尝试在代码中添加注释以使其更易于理解。看看 this dplyr cheat sheet.
基本上，您使用 group_by 根据 sample.number 将数据框分组。然后使用 summarise 计算每个组内 x 列的摘要指标。
要计算密度，您可以使用 summarise 中基数 R 的 density()。这将 return 一个包含密度函数 (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)
我有一个非常大的数据文件（>300k 行），每一行都是一个独特样本（>3000 个样本）的一部分。我想为每个单独的样本生成一个核密度估计器，并将相关信息（最小值、最大值、密度估计器的最大概率、密度估计器的中值、密度估计器的平均值）提取到单独的 table 以及样本名称。
我已尝试使用此处列出的方法从 ggplot 函数 stat_density_ridges() 中提取信息
Adding a mean to geom_density_ridges and here purrr::pluck 从 stat_density_ridges 和 ggplot_build 中提取数据，但它没有提供我想要的所有信息。
下面生成一些类似于我想要的合成数据：
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()
我想得到一个包含以下信息的数据 table：
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 中的管道，因此如果解决方案中包含管道，则可能需要一些简单的解释。
有多种方法可以做到这一点。在我看来，使用 dplyr 和管道运算符是最简单的方法。我尝试在代码中添加注释以使其更易于理解。看看 this dplyr cheat sheet.
基本上，您使用 group_by 根据 sample.number 将数据框分组。然后使用 summarise 计算每个组内 x 列的摘要指标。
要计算密度，您可以使用 summarise 中基数 R 的 density()。这将 return 一个包含密度函数 (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)