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

Question

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

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