R中具有两项（y=ax^b+cx^d）的幂函数的线性回归模型答案

【问题标题】：Linear regression model of a power function with two terms (y=ax^b+cx^d) in RR中具有两项（y=ax^b+cx^d）的幂函数的线性回归模型
【发布时间】：2020-07-09 17:25:57
【问题描述】：

我试图找出持续风速和阵风乘数之间的关系 [阵风乘数 =（阵风值）/（持续风速）]。我最初使用幂模型 (y = ax^b) 和线性回归对它们进行建模，但我想尝试将数据拟合到具有两个项 (y = ax^b + cx^d) 的幂模型中.

在具有一项的幂模型中，我能够获取每个变量的对数，然后我可以保存幂模型的系数（a 和 b 值）。我不确定如何更改我的代码以创建具有两个项的幂模型。

model = lm(log(data$GustMult)~log(data$SusWindSpeed))  # not sure how to change for a power model with two terms
coefs = unname(coef(model))
coefa = exp(coefs[1])
coefb = coefs[2]
coefc = exp(coefs[3])  #added for new model
coefd = coefs[4]  #added for new model
rsquared = summary(model)$r.squared

我尝试按照上面的方式运行代码，但它仍然只提供 a 和 b 系数的值。

提前感谢您提出的任何建议！

【问题讨论】：

欢迎堆栈溢出。如果您使您的问题可重现，则更容易提供帮助：以对象的形式包含一个最小数据集，例如，如果数据框为 df minimal reproducible example 和 How to Ask
我没有证据给你，但我强烈怀疑这没有被识别，这意味着你找不到唯一组参数（a，b， c，d）最适合该模型的数据。典型示例是尝试从 y = c + (a+b)x 中识别 a 和 b。在此示例中，存在无限数量的 a 和 b 组合，例如，可以最小化平方误差之和。我怀疑你提出的模型也有类似的问题。

标签： r linear-regression

【解决方案1】：

我认为 DanY 可能是对的，没有适合此模型的唯一参数集。无论如何，我认为您需要将其拟合为非线性模型。

我不知道你的数据是什么样的，但假设它是这样的：

set.seed(69)

SusWindSpeed <- rgamma(1000, 2, .125)
GustMult  <- (-0.98)*SusWindSpeed^(0.12) + (0.1)*SusWindSpeed^(-0.32) + rnorm(1000, 3, 0.025)
data <- data.frame(SusWindSpeed, GustMult)

plot(data$SusWindSpeed, data$GustMult)

然后您可以尝试像这样拟合非线性模型：

modA <- nls(GustMult ~ var_a * SusWindSpeed^var_b + var_c * SusWindSpeed^var_d + Const,
            start = list(var_a = -1, var_b = 0.1, var_c = 0.1, var_d = -0.3, Const = 3),
            data = data,
            control = list(maxiter = 500))

summary(modA)
#> 
#> Formula: GustMult ~ var_a * SusWindSpeed^var_b + var_c * SusWindSpeed^var_d + 
#>     Const
#> 
#> Parameters:
#>       Estimate Std. Error t value Pr(>|t|)    
#> var_a -1.28314    0.72635  -1.767   0.0776 .  
#> var_b  0.10381    0.04387   2.367   0.0181 *  
#> var_c  0.01506    0.03884   0.388   0.6983    
#> var_d -1.45958    3.32064  -0.440   0.6604    
#> Const  3.38496    0.76287   4.437 1.01e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.02452 on 995 degrees of freedom
#> 
#> Number of iterations to convergence: 22 
#> Achieved convergence tolerance: 3.379e-07

我们可以证明这是一个合理的拟合：

test_SusWindSpeed <- seq(0, 60, 0.05)

plot.new()
plot(data$SusWindSpeed, data$GustMult)
points(test_SusWindSpeed, 
       predict(modA, newdata = list(SusWindSpeed = test_SusWindSpeed)), 
       col = "red")

但是，您可能会发现您很难让该模型收敛，或者使用不同的起始参数会得到不同的结果。

^{由reprex package (v0.3.0) 于 2020-07-09 创建}

【讨论】：