Создание большого количества пользовательских функций в ggplot в R с использованием stat_function()

Question

Основная проблема заключается в том, что я хотел бы выяснить, как добавить большое количество (1000) пользовательских функций в один и тот же рисунок в ggplot, используя разные значения для коэффициентов функции. Я видел другие вопросы о том, как добавить две или три функции, но не 1000, и вопросы о добавлении в разных функциональных формах, но не в одной и той же форме с несколькими значениями параметров...

Цель состоит в том, чтобы stat_function рисовал линии, используя значения параметров, хранящиеся во фрейме данных, но без фактических данных для x.

[Общая цель здесь - показать большую неопределенность в параметрах модели нелинейной регрессии из небольшого набора данных, что выражается в неопределенности, связанной с предсказаниями на основе этих данных (что я пытаюсь убедить кого-то другого - плохая идея). Я часто делаю это, рисуя множество линий, построенных на основе неопределенности в параметрах модели (а-ля Эндрю Гелман, учебник по многоуровневой регрессии).]

В качестве примера приведем сюжет в базовой R графике.

#The data
p.gap <- c(50,45,57,43,32,30,14,36,51)
p.ag <- c(43,24,52,46,28,17,7,18,29)
data <- as.data.frame(cbind(p.ag, p.gap))

#The model (using non-linear least squares regression):
fit.1.nls <- nls(formula=p.gap~beta1*p.ag^(beta2), start=list(beta1=5.065, beta2=0.6168))
summary(fit.1.nls)

#From the summary, I find the means and s.e's the two parameters, and develop their distributions:
beta1 <- rnorm(1000, 7.8945, 3.5689)
beta2 <- rnorm(1000, 0.4894, 0.1282)
coefs <- as.data.frame(cbind(beta1,beta2))

#This is the plot I want (using curve() and base R graphics):
plot(data$p.ag, data$p.gap, xlab="% agricultural land use",
     ylab="% of riparian buffer gap", xlim=c(0,130), ylim=c(0,130), pch=20, type="n")
for (i in 1:1000){curve(coefs[i,1]*x^(coefs[i,2]), add=T, col="grey")}
curve(coef(fit.1.nls)[[1]]*x^(coef(fit.1.nls)[[2]]), add=T, col="red")
points(data$p.ag, data$p.gap, pch=20)

Я могу построить среднюю модель с данными в ggplot:

fit.mean <- function(x){7.8945*x^(0.4894)}
ggplot(data, aes(x=p.ag, y=p.gap)) +
  scale_x_continuous(limits=c(0,100), "% ag land use") +
  scale_y_continuous(limits=c(0,100), "% riparian buffer gap") +
  stat_function(fun=fit.mean, color="red") +
  geom_point()

Но ничто, что я делаю, не рисует несколько линий в ggplot. Кажется, я не могу найти какую-либо помощь в получении значений параметров из функций на веб-сайте ggplot или на этом сайте, которые обычно очень полезны. Это нарушает достаточно теории заговора, что никто не смеет делать это?

Любая помощь приветствуется. Спасибо!

Answer 1

Можно собрать несколько геомов или характеристик (и даже других элементов графика) в вектор или список и добавить этот вектор / список на график. Используя это, plyr пакет может быть использован для составления списка stat_functionпо одному на каждый ряд coefs

library("plyr")
coeflines <-
alply(as.matrix(coefs), 1, function(coef) {
  stat_function(fun=function(x){coef[1]*x^coef[2]}, colour="grey")
})

Затем просто добавьте это к сюжету

ggplot(data, aes(x=p.ag, y=p.gap)) +
  scale_x_continuous(limits=c(0,100), "% ag land use") +
  scale_y_continuous(limits=c(0,100), "% riparian buffer gap") +
  coeflines +
  stat_function(fun=fit.mean, color="red") +
  geom_point()

Пара заметок:

• Это медленно. На моем компьютере потребовалось несколько минут, чтобы нарисовать. ggplot не был разработан, чтобы быть очень эффективным при обработке около 1000 слоев.
• Это просто адрес добавления 1000 строк. За комментарий @ Роланда, я не знаю, представляет ли это то, что вы хотите / ожидаете статистически.

Answer 2

Вы могли бы создать новый stat_functions / изменить stat_function принять fun как эстетическое, как это:

# based on code from hadley and others
# found on https://github.com/tidyverse/ggplot2/blob/master/R/stat-function.r
library(rlang)
StatFunctions <- ggproto("StatFunctions", Stat,
                         default_aes = aes(y = stat(y)),
                         required_aes = "fun",

                         compute_group = function(data, scales, xlim = NULL, n = 101, args = list()) {
                           range <- xlim %||% scales$x$dimension()
                           xseq <- seq(range[1], range[2], length.out = n)

                           if (scales$x$is_discrete()) {
                             x_trans <- xseq
                           } else {
                             # For continuous scales, need to back transform from transformed range
                             # to original values
                             x_trans <- scales$x$trans$inverse(xseq)
                           }
                          do.call(rbind,
                                  lapply(data$fun, function(fun)
                                    data.frame(
                                      x = xseq,
                                      y =  do.call(fun, c(list(quote(x_trans)), args))))
                          )
                          }
)

stat_functions <- function(mapping = NULL, data = NULL,
                           geom = "path", position = "identity",
                           ...,
                           xlim = NULL,
                           n = 101,
                           args = list(),
                           na.rm = FALSE,
                           show.legend = NA,
                           inherit.aes = TRUE) {
  layer(
    data = data,
    mapping = mapping,
    stat = StatFunctions,
    geom = geom,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      n = n,
      args = args,
      na.rm = na.rm,
      xlim = xlim,
      ...
    )
  )
}

И затем используйте это так:

df <- data.frame(fun=1:3)
df$fun = c(function(x) x, function(x) x^2, function(x) x^3)
ggplot(df,aes(fun=fun, color=as.character(fun)))+
  stat_functions() +
  xlim(c(-5,5))

Чтобы получить это: