5 способов создания 2D гистограмм в R

==**Введение**==

Данная заметка является переводом [[http://www.everydayanalytics.ca/2014/09/5-ways-to-do-2d-histograms-in-r.html|http://www.everydayanalytics.ca/2014/09/5-ways-to-do-2d-histograms-in-r.html]].

[[http://www.mylesharrison.com/|Myles Harrison]] решил собрать как можно больше примеров создания качественных 2D гистограмм в R, для чего изучал различные форумы, блоги и задавал вопросы на [[http://stackoverflow.com/questions/18089752/r-generate-2d-histogram-from-raw-data|Stackoverflow]]. Эти примеры изложены ниже, с тем, чтобы их можно было использовать в работе, а также была наглядно видна разница между ними.

Для сведения, 2D гистограммы – это расширение обычных [[http://en.wikipedia.org/wiki/Histogram|гистограмм]], показывающих распределение частотности двух количественных показателей по некоторым интервалам группировки показателей. Их можно рассматривать как особый случай [[http://en.wikipedia.org/wiki/Heat_map|тепловых карт]], где цвет или интенсивность цвета показывает количество наблюдений попавших в соответствующий сегмент данных.

Во-первых, на протяжении всего текста будет использоваться цветовая палитра, созданная с помощью пакета [[http://cran.r-project.org/web/packages/RColorBrewer/index.html|RColorBrewer]], и нормально распределённые данные, объединённые в датафрейм. Во-вторых, для демонстрации полезности 2D гистограмм продемонстрируем соответствующую диаграмму рассеяния.

{{{ lang=rsplus
# Color housekeeping
library(RColorBrewer)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral'))) r <- rf(32) # Create normally distributed data for plotting x <- rnorm(mean=1.5, 5000) y <- rnorm(mean=1.6, 5000) df <- data.frame(x,y) # Plot plot(df, pch=16, col='black', cex=0.5) }}} [[image:scatterplot.jpeg|medium|center|link=source]] ==**Вариант 1: hexbin**== Пакет [[http://cran.r-project.org/web/packages/hexbin/index.html|hexbin]] разбивает пространство на 2D шестигранные соты и затем подсчитывает количество попаданий в каждую соту. Удобством использования пакета //hexbin// является то, что он автоматически добавляет легенду в график, так как ручное добавление легенды в R вещь достаточно проблематичная. Если использовать функции с параметрами по умолчанию, то получится простой монохромный график. Использование параметра //colramp// с нужной палитрой цветов немного улучшает картинку. К сожалению, положение легенды оставляет желать лучшего. {{{ lang=rsplus #----- OPTION 1: hexbin from package 'hexbin' -----------[[image:kde2d_1.jpeg|center|large|link=source]] library(hexbin) # Create hexbin object and plot h <- hexbin(df) plot(h) plot(h, colramp=rf) }}} [[image:hexbin_1.jpeg|large|center|link=source]] Функции //hexbinplot// дает большую гибкость, позволяя настройку границ бинирования, к тому же поддерживает использование функций трансформации, к примеру, на одном из нижеследующих графиков была использована логарифмическая шкала. Только не стоит забывать, что необходимо передать как прямую, так и обратную функцию трансформации. Дополнительным плюсом данной функции является то, что корректно вычисляется положение легенды. {{{ lang=rsplus # hexbinplot function allows greater flexibility hexbinplot(y~x, data=df, colramp=rf) # Setting max and mins hexbinplot(y~x, data=df, colramp=rf, mincnt=2, maxcnt=60) }}} [[image:hexbin_2.jpeg|large|center|link=source]] {{{ lang=rsplus # Scaling of legend - must provide both trans and inv functions hexbinplot(y~x, data=df, colramp=rf, trans=log, inv=exp) }}} [[image:hexbin_3-e1409740343841.jpeg|center|link=source]] ==**Вариант 2: hist2d**== Пакет [[http://cran.r-project.org/web/packages/gplots/index.html|gplots]] помимо всего прочего содержит функцию //hist2d//, которая позволяет быстро строить 2D гистограммы. И вновь, функция с параметрами по умолчанию даёт не очень эффектный график {{{ lang=rsplus #----- OPTION 2: hist2d from package 'gplots' --------- library(gplots) # Default call h2 <- hist2d(df) }}} [[image:hist2d_1.jpeg|center|medium|link=source]] Настройка цветовой палитры и размера сегментов даёт более хороший результат. И так же, как и в предыдущем случае можно настроить логарифмическую шкалу. {{{ lang=rsplus # Coarser binsizing and add colouring h2 <- hist2d(df, nbins=25, col=r) # Scaling with log as before h2 <- hist2d(df, nbins=25, col=r, FUN=function(x) log(length(x))) }}} [[image:hist2d_adj.jpeg|center|large|link=source]] ==**Вариант 3: stat_2bin из пакета ggplot2**== Конечно же, разве может какая-либо статья по R обойтись без упоминания пакета [[http://cran.r-project.org/web/packages/ggplot2/index.html|ggplot2]]? При построении 2D гистограмм в этом пакете можно использовать функцию //stat_2bin// либо добавляя её к //ggplot//-объекту, либо передавая в качестве параметра типа геометрии в функции //qplot//. {{{ lang=rsplus #----- OPTION 3: stat_bin2d from package 'ggplot' ------- library(ggplot2) # Default call (as object) p <- ggplot(df, aes(x,y)) (h3 <- p + stat_bin2d()) # Default call (using qplot) qplot(x,y,data=df, geom='bin2d') }}} [[image:ggplot_1-e1409741740122.jpeg|center|link=source]] И вновь, имеет смысл поменять сегментирование и использовать свою цветовую палитру. Последнее можно сделать при помощи функции //scale_fill_gradientn// с нужным вектором цветов. Логарифмическую шкалу можно задать в той же функции, передав её в параметр //trans//. {{{ lang=rsplus # Add colouring and change bins (h3 <- p + stat_bin2d(bins=25) + scale_fill_gradientn(colours=r)) # Log scaling (h3 <- p + stat_bin2d(bins=25) + scale_fill_gradientn(colours=r, trans="log")) }}} [[image:ggplot_2.jpeg|center|large|link=source]] ==**Вариант 4: kde2d**== В этом варианте используется ядерная оценка плотности (англ. kernel denisty estimation) при помощи функции //kde2d// из пакета [[http://cran.r-project.org/web/packages/MASS/index.html|MASS]]. Здесь стоит учесть, что получаемый результат не является в точности 2D гистограммой, так как функция интерполирует данные для получения соответствующей оценки плотности распределения. По умолчанию используется значение количества сегментов n = 25, что совпадает с тем, что использовалось ранее, а потому этот параметр не будет указываться явно. Изображение строится при помощи функции //image()//. Увеличение числа n фактически будет приближать результат к ядерной оценке плотности, так как в какой-то момент размер сегмента станет меньше, чем шаг в имеющихся данных. Hadley Wickham написал статью, в которой указал, что в R [[http://vita.had.co.nz/papers/density-estimation.pdf|более 20 пакетов (PDF)]], при помощи которых можно строить оценки плотности распределения, поэтому мы не будем углубляться в эту тему. {{{ lang=rsplus #----- OPTION 4: kde2d from package 'MASS' ------------- # Not a true heatmap as interpolated (kernel density estimation) library(MASS) # Default call k <- kde2d(df$x, df$y) image(k, col=r) # Adjust binning (interpolate - can be computationally intensive for large datasets) k <- kde2d(df$x, df$y, n=200) image(k, col=r) }}} [[image:kde2d_1.jpeg|center|large|link=source]] ==**Вариант 5: Использование базовых библиотек**== И, наконец, бесстрашный пользователь R продемонстрировал на [[http://stackoverflow.com/questions/18089752/r-generate-2d-histogram-from-raw-data/18103689#18103689|Stackoverflow]] как можно решить эту задачу "в лоб", с использованием только базовых библиотек. {{{ lang=rsplus #----- OPTION 5: The Hard Way (DIY) ----------------- # http://stackoverflow.com/questions/18089752/r-generate-2d-histogram-from-raw-data nbins <- 25 x.bin <- seq(floor(min(df[,1])), ceiling(max(df[,1])), length=nbins) y.bin <- seq(floor(min(df[,2])), ceiling(max(df[,2])), length=nbins) freq <- as.data.frame(table(findInterval(df[,1], x.bin),findInterval(df[,2], y.bin))) freq[,1] <- as.numeric(freq[,1]) freq[,2] <- as.numeric(freq[,2]) freq2D <- diag(nbins)*0 freq2D[cbind(freq[,1], freq[,2])] <- freq[,3] # Normal image(x.bin, y.bin, freq2D, col=r) # Log image(x.bin, y.bin, log(freq2D), col=r) }}} [[image:hardway_1.jpeg|center|large|link=source]] Не тот способ, которым стоит это делать, особенно с учётом вышеописанных возможностей, но если есть настоятельная необходимость не использовать дополнительные пакеты, то это рабочий вариант. ==**Бонусный вариант**== И напоследок, есть смысл добавить эту, очень хорошую диаграмму из книги [[http://books.google.ca/books?id=YWcLBAAAQBAJ&printsec=frontcover&dq=Computational+Actuarial+Science+with+R&hl=en&sa=X&ei=0hYHVOhrgZrKA5S9gLAI&ved=0CDEQ6AEwAA#v=onepage&q=Computational%20Actuarial%20Science%20with%20R&f=false|Computational Actuarial Science with R]]. Эта диаграмма используется не очень часто и состоит из 2D гистограммы и двух обычных 1D гистограмм, показывающих распределение величины в каждом из измерений. {{{ lang=rsplus #----- Addendum: 2D Histogram + 1D on sides (from Computational ActSci w R) ------ #http://books.google.ca/books?id=YWcLBAAAQBAJ&pg=PA60&lpg=PA60&dq=kde2d+log&source=bl&ots=7AB-RAoMqY&sig=gFaHSoQCoGMXrR9BTaLOdCs198U&hl=en&sa=X&ei=8mQDVPqtMsi4ggSRnILQDw&redir_esc=y#v=onepage&q=kde2d%20log&f=false h1 <- hist(df$x, breaks=25, plot=F) h2 <- hist(df$y, breaks=25, plot=F) top <- max(h1$counts, h2$counts) k <- kde2d(df$x, df$y, n=25) # margins oldpar <- par() par(mar=c(3,3,1,1)) layout(matrix(c(2,0,1,3),2,2,byrow=T),c(3,1), c(1,3)) image(k, col=r) #plot the image par(mar=c(0,2,1,0)) barplot(h1$counts, axes=F, ylim=c(0, top), space=0, col='red') par(mar=c(2,0,0.5,1)) barplot(h2$counts, axes=F, xlim=c(0, top), space=0, col='red', horiz=T) }}} [[image:addendum.jpeg|center|medium|link=source]] ==**Заключение**== Итак, было предоставлено 5 различных способов построения 2D гистограмм в R, плюс некоторый дополнительный код. В качестве самостоятельного задания остаётся задача добавления легенд в те гистограммы, в которых её нет. Надеюсь, что пользователи R найдут эту заметку полезной. Весь код, описанный в тексте можно скачать [wpfilebase tag="fileurl" id=10 linktext="здесь"/].