【问题标题】:Power regression in R similar to excelR中的幂回归类似于excel
【发布时间】:2018-05-08 09:22:46
【问题描述】:

我有一个简单的数据集,我正在尝试使用功率趋势来最好地拟合数据。样本数据很小,如下:

structure(list(Discharge = c(250, 300, 500, 700, 900), Downstream = c(0.3, 
0.3, 0.3, 0.3, 0.3), Age = c(1.32026239202165, 1.08595138888889, 
0.638899189814815, 0.455364583333333, 0.355935185185185)), .Names = c("Discharge", 
"Downstream", "Age"), row.names = c(NA, 5L), class = "data.frame")

数据如下:

> new
  Discharge Downstream       Age
1       250        0.3 1.3202624
2       300        0.3 1.0859514
3       500        0.3 0.6388992
4       700        0.3 0.4553646
5       900        0.3 0.3559352

我尝试使用ggplot2 绘制上述数据

ggplot(new)+geom_point(aes(x=Discharge,y=Age))

我可以使用geom_smooth(method="lm") 添加线性线,但我不确定需要什么代码来显示电源线。

输出如下:

如何像在 excel 中那样添加幂线性回归线? excel图如下:

【问题讨论】:

    标签: r ggplot2


    【解决方案1】:

    虽然 mnel 的答案对于非线性最小二乘拟合是正确的,但请注意,Excel 实际上并没有做任何几乎那么复杂的事情。它实际上只是对响应变量和预测变量进行对数转换,并进行普通(线性)最小二乘拟合。要在 R 中重现这一点,您可以:

    lm(log(Age) ~ log(Discharge), data=df)
    
    Call:
    lm(formula = log(Age) ~ log(Discharge), data = df)
    
    Coefficients:
       (Intercept)  log(Discharge)  
             5.927          -1.024  
    

    作为检查,log(Discharge) 的系数与 Excel 中的系数相同,而 exp(5.927) ~ 375.05。

    虽然我不确定如何在 ggplot2 中将其用作趋势线,但您可以在基本图形中这样做:

    m <- lm(log(y) ~ log(x), data=df)
    
    newdf <- data.frame(Discharge=seq(min(df$Discharge), max(df$Discharge), len=100))
    plot(Age ~ Discharge, data=df)
    lines(newdf$Discharge, exp(predict(m, newdf)))
    
    text(600, .8, substitute(b0*x^b1, list(b0=exp(coef(m)[1]), b1=coef(m)[2])))
    text(600, .75, substitute(plain("R-square: ") * r2, list(r2=summary(m)$r.squared)))
    

    【讨论】:

    • 确实——因此我对发布 r^2 值保持沉默——
    • 使用ggplot2 这可以使用stat_smooth(method = 'glm', formula = 'y~log(x)', family = gaussian(link = 'log') 进行拟合
    • @mnel 但这并不完全相同。这适合log(E(Y)) = log(x),而不是E(log(Y)) = log(x)。我不认为有一种方法可以将对数变换回归拟合为 ggplot2 中的平滑 - 这是有道理的,因为回归在原始尺度上存在偏差。
    • 好点。也许scale_x_continuous(trans = 'log') + scale_y_continuous(trans = 'log') + stat_smooth(method = 'lm', se = FALSE) + coord_trans(ytrans = 'exp', xtrans = 'exp')(这会记录数据并将平滑(线性)拟合到转换后的数据,然后反转转换以进行绘图。
    • 对当前版本的 ggplot 的 @mnel 答案进行了细微更改。它应该生成匹配excel的ggplot曲线:ggplot(df,aes(x,y)) + geom_point() + scale_x_continuous(trans = 'log') + scale_y_continuous(trans = 'log') + stat_smooth(method = 'lm', se = FALSE) + coord_trans(y = 'exp', x = 'exp')
    【解决方案2】:

    使用nls(非线性最小二乘法)作为平滑器

    例如

    ggplot(DD,aes(x = Discharge,y = Age)) +
      geom_point() + 
      stat_smooth(method = 'nls', formula = 'y~a*x^b', start = list(a = 1,b=1),se=FALSE)
    

    注意 Doug Bates cmets 关于 R 平方值和非线性模型 here,您可以使用以下想法 Adding Regression Line Equation and R2 on graph

    附加回归线方程

    # note that you have to give it sensible starting values
    # and I haven't worked out why the values passed to geom_smooth work!
    power_eqn = function(df, start = list(a =300,b=1)){
      m = nls(Discharge ~ a*Age^b, start = start, data = df);
      eq <- substitute(italic(y) == a  ~italic(x)^b, 
                   list(a = format(coef(m)[1], digits = 2), 
                        b = format(coef(m)[2], digits = 2)))
      as.character(as.expression(eq));                 
    }
    
    ggplot(DD,aes(x = Discharge,y = Age)) +
      geom_point() + 
      stat_smooth(method = 'nls', formula = 'y~a*x^b', start = list(a = 1,b=1),se=FALSE) +  
      geom_text(x = 600, y = 1, label = power_eqn(DD), parse = TRUE)
    

    【讨论】:

    • 非常感谢您的快速回答。您是否还建议我如何在ggplot 之外拟合非线性模型来发布R2 的值和图中的方程?
    • 我使用了以下公式nlm &lt;- nls(Age ~ a*Discharge^b,data=new,start=list(a=1,b=1))。你能让我知道它是否正确吗?我似乎无法使用公式找到R-squared
    • @Jdbaba -- 你需要合理的起始值(尝试 a = 300),我还没有弄清楚为什么 (a=1) 不会在 geom_smooth 中抛出相同的错误(很奇怪)。我添加了一个关于 R^2 和非线性最小二乘法的链接。
    【解决方案3】:

    2018 年更新: 电话"start" 现在似乎已经贬值了。 stat_smooth 函数信息中也没有。

    如果要选择起始值,现在需要使用“method.args”选项。

    查看以下更改:

    ggplot(DD,aes(x = Discharge,y = Age)) +
      geom_point() + 
      stat_smooth(method = 'nls', formula = 'y~a*x^b', method.args = list(start= c(a = 1,b=1)),se=FALSE) + geom_text(x = 600, y = 1, label = power_eqn(DD), parse = TRUE)
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2012-12-21
      • 1970-01-01
      • 1970-01-01
      • 2018-03-18
      • 2015-11-10
      • 2011-03-27
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多