Объедините игровые сглаживания нескольких кластеров на графике mgcViz

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

Пакет mgcViz обеспечивает отличную визуализацию для отдельных кластеров, но я не вижу, как их объединить. Возможно, это потому, что он направлен на визуализацию нескольких эффектов, а не нескольких кластеров. Тем не менее, его возможности очень близки к тому, что мне нужно, поэтому вот воспроизводимый пример (адаптировано из https://mfasiolo.github.io/mgcViz/articles/mgcviz.html):

library(mgcViz)
n = 1e3
z = rnorm(n)
dat = data.frame(x = rep(z, times = 2),
                 y = rep(c(1,2), each = n) + c(sin(z), 0.5*z^2) + rnorm(2*n)/4,
                 g = factor(rep(1:2, each = n)))

b <- lapply(1:2, function(i, dat) gam(y ~ s(x), data = dat[dat$g == i, ]),
            dat = dat)

plot(getViz(b[[1]])) + l_points() + l_fitLine() + l_ciLine()   # First
plot(getViz(b[[2]])) + l_points() + l_fitLine() + l_ciLine()   # Second
plot(getViz(b))   # Third
ggplot(dat, aes(x, y, color = g)) + geom_point(pch = ".") + theme_bw() # Fourth

Первый сюжетВторой сюжет Третий сюжетЧетвертый сюжет

Я хотел бы объединить первые два сюжета в один, как это частично сделано в третьем сюжете. Включение третьего графика в данные, показанные в четвертом, поможет. Это также требует различных сдвигов перехвата в третьем графике. Добавление l_points() к третьему графику делает его пустым.

Скрытое ограничение заключается в том, что сглаживания gam являются отдельными компонентами списка (как показано выше), потому что они фактически исходят из пользовательского кода кластеризации для фрагментов временных рядов разной длины и интервалов с использованием bam mgcv для очень больших данных. Таким образом, график предпочтительно должен получать всю информацию из b, список результатов gam для каждого кластера.


r cluster-analysis ggplot2 mgcv
person George Ostrouchov    schedule 28.01.2021    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Не mgcViz, но вы можете просто создать необходимый вывод самостоятельно и использовать ggplot2:

library(mgcv)
library(ggplot2)
theme_set(theme_bw())

n = 1e3
z = rnorm(n)
dat = data.frame(
  z = rep(z, times = 2),
  y = rep(c(1,2), each = n) + c(sin(z), 0.5*z^2) + rnorm(2*n)/4,
  g = factor(rep(1:2, each = n)))

b <- gam(y ~ g + s(z, by = g), data = dat)

ndf <- expand.grid(z = seq(min(dat$z), max(dat$z), length.out=100), g = unique(dat$g))
ndf$pred <- predict(b, newdata = ndf, type = "response")

ggplot(ndf, aes(x = z, y = pred, col = g)) +
  geom_line() +
  geom_point(data = dat, aes(y = y))

Создано 28 января 2021 г. с помощью пакета reprex (v0.3.0)

Редактировать

Если вы хотите поместить отдельные модели в группу, можно сделать что-то вроде этого:

library(purrr)
ndf <- map_dfr(
  .x = unique(dat$g),
  .f = ~{
    mod_i <- gam(y ~ s(x), data = dat[dat$g == .x, ])
    ndf_i <- expand.grid(x = seq(min(dat$x), max(dat$x), length.out = 100))
    ndf_i$g <- .x
    ndf_i$pred <- predict(mod_i, newdata = ndf_i, type = "response")
    ndf_i
  })

ggplot(ndf, aes(x = x, y = pred, col = g)) +
  geom_line() +
  geom_point(data = dat, aes(y = y))
person adibender    schedule 28.01.2021
comment
Мой b представляет собой список отдельных результатов gam для каждого кластера (я обновил вопрос), поэтому у меня нет ни одной модели с g. Я мог бы переделать модель с перехватами g, но мне интересно, можно ли этого избежать, поскольку задействованные данные велики. - person George Ostrouchov; 29.01.2021
comment
Это работает, с немного большим количеством данных, если нужно начать с b, поскольку список результатов gam в моем вопросе. В частности, ваш predict заменен на do.call(c, lapply(b, predict, newdata = list(x = z), type = "response")). Спасибо! - person George Ostrouchov; 29.01.2021
comment
@GeorgeOstrouchov Я внес изменения в сообщение по запросу. Обратите внимание, что в версии с одной моделью вы также можете использовать функцию mgcv::bam вместо gam с method = "fREML" и discrete = TRUE, чтобы сократить время вычислений и требования к памяти. - person adibender; 29.01.2021