Как получить уравнение для линии регрессии на логарифмическом графике в ggplot2?

У меня есть логарифмический график, я получил линию регрессии, используя:

geom_smooth(formula = y ~ x, method='lm') 

Но теперь я хотел бы получить уравнение этой линии (например, y=a*x^(-b)) и распечатать его. Мне удалось получить это в сюжете лин-лин, но не в этом случае. Вот код:

mydataS<-data.frame(DurPeak_h[],IntPeak[],IntPeakxDurPeak[],ID[]) #df peak
names(mydataS)<-c("x","y","ID","IDEVENT")

plotID<-ggplot(mydataS, aes(x=x, y=y, label=IDEVENT)) + 
geom_text(check_overlap = TRUE, hjust = 0, nudge_x = 0.02)+
geom_point(colour="black", size = 2) + geom_point(aes(colour = ID)) +
geom_quantile(quantiles = qs, colour="green")+ 
scale_colour_gradient(low = "white", high="red") +
scale_x_log10(limits = c(min(DurEnd_h),max(DurEnd_h))) + 
scale_y_log10(limits = c(min(IntEnd),max(IntEnd))) +
geom_smooth(formula = y ~ x, method='lm') 

ggsave(height=7,"plot.pdf")

person chil989    schedule 04.01.2017    source источник
comment
вероятно, проще всего подогнать модель вне ggplot   -  person Richard Telford    schedule 04.01.2017
comment
Как вы думаете, откуда взялось это уравнение?   -  person Roman Luštrik    schedule 04.01.2017
comment
Не усложняйте себе жизнь больше, чем она должна быть, просто подгоните lm под себя. ggplot2 не является пакетом моделирования.   -  person Axeman    schedule 04.01.2017


Ответы (2)


mydataS<-data.frame(DurPeak_h[],IntPeak[],IntPeakxDurPeak[],ID[])
names(mydataS)<-c("x","y","ID","IDEVENT")
model <- lm(y~x, header = T)
summary(model)

используйте значение перехвата, указанное как «b», и коэффициент как ваш «a»

person Amy    schedule 04.01.2017

Сделал это обходным путем: используя nls для вычисления двух параметров a и b, а именно:

nlsPeak <- coef(nls(y ~ a*(x)^b, data = mydataS, start = list(a=30, b=-0.1)))

затем построить линию с annotate (см. некоторые примеры здесь) и, наконец, распечатать уравнение, используя функция:

power_eqn = function(ds){
  m = nls(y ~ a*x^b, start = list(a=30, b=-0.1), data = ds);
  eq <- substitute(italic(y) == a  ~italic(x)^b, 
               list(a = format(coef(m)[1], digits = 4), 
                    b = format(coef(m)[2], digits = 2)))
  as.character(as.expression(eq));
}

называется следующим образом:

annotate("text",x = 3, y = 180,label = power_eqn(mydataS), parse=TRUE, col="black") +

Надеюсь, поможет!

person chil989    schedule 05.01.2017