Как я могу извлечь данные из функции плотности ядра в R для многих образцов за один раз

Question

У меня очень большой файл данных (> 300 тыс. Строк), и каждая строка является частью уникальной выборки (> 3000 выборок). Я хочу сгенерировать оценку плотности ядра для каждого отдельного образца и извлечь соответствующую информацию (минимальное значение, максимальное значение, максимальная вероятность оценки плотности, медиана оценки плотности, среднее значение оценки плотности) в отдельную таблицу вместе с названием образца.

Я попытался извлечь информацию из функции ggplot stat_density_ridges(), используя изложенные здесь подходы Добавление среднего значения для geom_density_ridges и здесь нарисовать линию на 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 и pipe-оператора, на мой взгляд, самый простой способ. Я попытался добавить комментарии в код, чтобы было легче понять. Взгляните на этот шпаргалки dplyr .

По сути, вы используете group_by, чтобы разделить ваш фрейм данных на группы в соответствии с sample.number. Затем вы используете summarise для вычисления итоговых метрик столбца x внутри каждой группы.

Для вычисления плотности вы можете использовать density() из базы R внутри 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)