Методы создания тепловых карт, графиков уровня / контура и гексагонального биннинга

Question

Вариантов для 2D-графиков (x, y, z) в R немного. Тем не менее, борьба с опциями - это сложная задача, особенно в том случае, если все три являются непрерывными.

Чтобы прояснить проблему (и, возможно, помочь объяснить, почему меня могут споткнуть с contour или image), вот возможная схема классификации:

• Случай 1: значение z не указано, но является условной плотностью, основанной на значениях в (x, y). (Примечание: это по существу сводит вычисление z к отдельной функции - оценке плотности. Что-то еще должно использовать результаты этого вычисления, поэтому было бы неплохо учесть произвольные вычисления.)
• Случай 2: (x, y) пары уникальны и расположены через равные интервалы. Это означает, что для каждого значения (x, y) предоставляется только одно значение z.
• Случай 3: (x, y) пары уникальны, но непрерывны. Цвет или затенение по-прежнему определяется только 1 уникальным значением z.
• Случай 4: (x, y) пары не уникальны, а располагаются регулярно. Цвет или затенение определяется функцией агрегирования по значениям z.
• Случай 5: (x, y) пары не являются уникальными, являются непрерывными. Цвет / затенение должны определяться функцией агрегирования по значениям z.

Если я пропускаю некоторые случаи, пожалуйста, дайте мне знать. Случай, который меня интересует, # 5. Некоторые заметки об отношениях:

• Дело № 1 уже хорошо поддерживается.
• Дело № 2 легко поддерживается heatmap, image и функционирует в ggplot.
• Случай № 3 поддерживается базой plot, хотя использование цветового градиента оставлено на усмотрение пользователя.
• Дело № 4 может стать делом № 2 с помощью функции разделения и применения. Я делал это раньше.
• Дело № 5 можно преобразовать в # 4 (а затем # 2) с помощью cut, но это не элегантно и квадратно. Шестнадцатеричное биннинг может быть лучше, хотя, кажется, это не может быть легко обусловлено наличием крутого градиента в значении z. Я бы согласился на шестнадцатеричное бинирование, но альтернативные функции агрегирования вполне приветствуются, особенно если они могут использовать значения z.

Как я могу сделать # 5? Вот код для создания седла, хотя значение spread изменяет разброс значения z, что должно создавать различия в построении градиентов.

N       = 1000
spread  = 0.6   # Vals: 0.6, 3.0
set.seed(0)
rot     = matrix(rnorm(4), ncol = 2)
mat0    = matrix(rnorm(2 * N), ncol = 2)
mat1    = mat0 %*% rot
zMean   = mat0[,2]^2 - mat0[,1]^2
z       = rnorm(N, mean = zMean, sd = spread * median(abs(zMean)))

Я бы хотел сделать что-то вроде hexbin, но я справился с этим с ggplot и не добился большого прогресса. Если бы я мог применить произвольную функцию агрегирования к значениям z в регионе, это было бы еще лучше. (Форма такой функции может быть такой, как plot(mat1, colorGradient = f(z), aggregation = "bin", bins = 50).)

Как я могу сделать это в ggplot или другом пакете? Я рад сделать этот вопрос вопросом вики сообщества (или другим пользователям, отредактировав его достаточно много раз). Если да, один ответ на сообщение, пожалуйста, чтобы мы могли сосредоточиться, скажем, на ggplot, levelplot, lattice, contourplot (или image) и других параметрах, если они существуют.

---

Обновления 1: Пример вулкана является хорошим примером случая № 3: данные регулярно располагаются (это может быть широта / долгота), с одним значением z на наблюдение. Топографическая карта имеет (широту, долготу, высоту) и, следовательно, одно значение на местоположение. Предположим, кто-то получает погоду (например, количество осадков, скорость ветра, солнечный свет) в течение многих дней для многих случайно расположенных датчиков: это больше похоже на № 5, чем на № 3 - у нас может быть lat & long, но значения z могут варьироваться совсем немного , даже для тех же или близких (x, y) значений.

Обновление 2: до сих пор ответы от DWin, Kohske и John Colby превосходны. Мой фактический набор данных представляет собой небольшую выборку из большего набора, но при 200 тыс. Точек он дает интересные результаты. На плоскости (x, y) он имеет очень высокую плотность в некоторых регионах (таким образом, в этих областях может произойти переполнение) и гораздо меньшую плотность или полное отсутствие в других регионах. С предложением Джона через fields мне нужно было отобрать данные для обработки Tps (я выясню, смогу ли я сделать это без подвыборки), но результаты довольно интересные. При попытке rms / Hmisc (предложение DWin) полные 200К точек, кажется, работают хорошо. Предложение Кохске довольно хорошее, и, поскольку данные преобразуются в сетку перед построением графика, нет проблем с количеством точек входных данных. Это также дает мне большую гибкость, чтобы определить, как агрегировать значения z в регионе. Я еще не уверен, буду ли я использовать среднее, медиану или какую-либо другую совокупность.

Я также собираюсь опробовать хороший пример Кохске с mutate + ddply с другими методами - это хороший пример того, как получить различную статистику, рассчитанную для данного региона.

---

Обновление 3: различные методы различны, а некоторые замечательны, хотя нет явного победителя. Я выбрал ответ Джона Колби первым. Я думаю, что я буду использовать этот метод или метод DWin в дальнейшей работе.

Answer 1

Мне очень повезло с пакетом fields для этого типа проблемы.Вот пример использования Tps для тонких сплайнов:

EDIT: объединенные графики и добавленная стандартная ошибка

require(fields)

dev.new(width=6, height=6)
set.panel(2,2)

# Plot x,y
plot(mat1)

# Model z = f(x,y) with splines
fit = Tps(mat1, z)
pred = predict.surface(fit)

# Plot fit
image(pred)
surface(pred)

# Plot standard error of fit 
xg = make.surface.grid(list(pred$x, pred$y))
pred.se = predict.se(fit, xg)

surface(as.surface(xg, pred.se))

Answer 2

В пакете latticeExtra есть функция panel.2dsmoother:

library(lattice)
library(latticeExtra)

df <- data.frame(mat1, z)
names(df)[c(1,2)] <- c('x', 'y')

levelplot(z ~ x * y, data = df, panel = panel.2dsmoother, contour=TRUE)

Согласно справочной странице"модель сглаживания строится (приблизительно) как метод (форма, данные = список (x = x, y = y, z = z), {args}) [ ...] Это должно работать с любой функцией модели, которая принимает аргумент формулы и имеет аргумент метода предиката ".

Answer 3

Я обычно использую комбинацию rms / Hmisc.Это линейный регрессионный анализ (функция ols) с использованием скрещенных кубических сплайн-членов, чей вывод на графике очень похож на предлагаемый пример полей:

dfrm <- data.frame(z=z, xcor = mat1[,1], ycor=mat1[,2])
require(rms)  # will automatically load Hmisc which needs to have been installed
lininterp <- ols(z ~ rcs(xcor,3)*rcs(ycor,3), data=dfrm)
ddI <- datadist(dfrm)
options(datadist="ddI")

 bplot(Predict(lininterp, xcor, ycor))  # Plot not shown
 perim <- with(dfrm, perimeter(xcor, ycor))
 bplot(Predict(lininterp, xcor, ycor), perim=perim)  
# Plot attached after converting to .png

Также можно увидеть метод, который не полагается на регрессионные оценки3-D поверхность во второй части моего ответа на этот вопрос: Использование цвета в качестве 3-го измерения

Парадигма построения решетки, и вы также можете получить контурные диаграммыкак этот симпатичный уровень.Если вы хотите, чтобы прогнозируемые значения находились в начальной точке, вы можете получить это с помощью функции Predict, примененной к объекту подгонки.

Answer 4

Вероятно, вопрос можно разделить на две части.Первый - это агрегирование данных, а второй - их визуализация. Пакет

fields, как показывает @John, может выполнять эти действия одновременно.В ggplot2, если агрегация - это просто подсчет данных, stat_bin2d доступно.

В любом случае, если вы хотите использовать собственную функцию агрегирования, может быть что-то вроде этого поможет:

df <- data.frame(x = mat1[,1], y = mat1[,2], z = z)

Nx <- 10 # nubmer of bins for x
Ny <- 4  # number of bins for y

# create a data.
df2 <- mutate(ddply(df, .(x = cut(x, Nx), y = cut(y, Ny)), summarise, 
                    Mean = mean(z),
                    Var = var(z)),
              xmin = as.numeric( sub("\\((.+),.*", "\\1", x)),
              xmax = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", x)),
              ymin = as.numeric( sub("\\((.+),.*", "\\1", y)),
              ymax = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", y)),
              xint = as.numeric(x),
              yint = as.numeric(y))

# then, visualize
ggplot(df2, aes(xint, yint, xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax, fill = Mean)) +
  geom_tile(stat = "identity")

ggplot(df2, aes(xint, yint, xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax, fill = Var)) +
  geom_tile(stat = "identity")