如何在 R 中拟合 a < 0 的二次模型？

Question

我正在对沿海拔梯度的蜜蜂多样性拟合二次模型。我假设沿梯度某处会有最大值，因此我的模型应该有一个负 "a" 系数。这适用于 3 个属，但对于第四个属 (Exaerete)，"a" 变为阳性。下图显示了所有 4 个拟合，我们可以看到蓝线是唯一的 "incorrect":

分离这个属，我们可以清楚地看到为什么它是"incorrect":

同时存在二次模型和线性模型。考虑到数据点，二次方有意义，但在生物学上意义不大。我想强制命令生成负数 "a"（因此 "optimum" 高度可能比第一张图中给定的高度低得多，即 1193 米），我该怎么做？ R 中用于生成模型的命令是

fitEx2 <- lm(num~I(alt^2)+alt,data=Ex)

而数据是

Ex <- data.frame(alt=c(50,52,100,125,130,200,450,500,525,800,890,1140),
                 num=c(3,1,2,1,1,2,1,2,1,1,1,1))

Answer 1

我们正在处理一个受限的估计，可以方便地处理，例如 nls。例如，

x <- rnorm(100)
y <- rnorm(100) - 0.01 * x^2 + 0.1 * x

nls(y ~ -exp(a) * x^2 + b * x + c, start = list(a = log(0.01), b = 0.1, c = 0))
# Nonlinear regression model
#   model: y ~ -exp(a) * x^2 + b * x + c
#    data: parent.frame()
#        a        b        c 
# -4.66893 -0.03615 -0.01949 
#  residual sum-of-squares: 97.09
# 
# Number of iterations to convergence: 2 
# Achieved convergence tolerance: 3.25e-08

其中使用 exp 有助于施加负约束。那么你想要的二次项系数是

-exp(-4.66893)
[1] -0.009382303

但是，由于 lm 估计的系数为正，因此在您的特定情况下，nls 可能会因接近 −∞ 以使系数为零而崩溃。

可能使用更稳定的方法，例如 optim:

set.seed(2)
x <- rnorm(100)
y <- rnorm(100) - 0.01 * x^2 + 0.1 * x
lm(y ~ x + I(x^2))

# Call:
# lm(formula = y ~ x + I(x^2))

# Coefficients:
# (Intercept)            x       I(x^2)  
#    -0.04359      0.04929      0.04343  

fun <- function(b) sum((y - b[1] * x^2 - b[2] * x - b[3])^2)
optim(c(-0.01, 0.1, 0), fun, method = "L-BFGS-B",
      lower = c(-Inf, -Inf, -Inf), upper = c(0, Inf, Inf))
# $par
# [1] 0.00000000 0.05222262 0.01441276
# 
# $value
# [1] 95.61239
# 
# $counts
# function gradient 
# 7        7 
# 
# $convergence
# [1] 0
# 
# $message
# [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

这表明线性模型。事实上，由于您的模型非常简单，很可能这确实是理论上的最优值，您可能需要重新考虑您的方法。例如，也许您可以将一些观察值视为异常值并相应地改变估计值？

Answer 2

Optimization Task View 列出了几个解决 least-squares 问题的包，允许对系数进行（线性）约束，例如 minpack.lm。

library(minpack.lm)
x <- Ex$alt; y <- Ex$num
nlsLM(y ~ a*x^2 + b*x + c, 
      lower=c(-1, 0, 0), upper=c(0, Inf, Inf), 
      start=list(a=-0.01, b=0.1, c=0))
## Nonlinear regression model
##   model: y ~ a * x^2 + b * x + c
##    data: parent.frame()
##     a     b     c 
## 0.000 0.000 1.522 
##  residual sum-of-squares: 5.051
## 
## Number of iterations to convergence: 27 
## Achieved convergence tolerance: 1.49e-08

顺便说一句，这个功能也比nls更可靠，并尽量避免"zero gradient"。
如果用户更频繁地利用许多 CRAN 任务视图，将会很有帮助。

如何在 R 中拟合 a < 0 的二次模型？

How to fit a quadratic model with a < 0 in R?

r

quadratic

curve-fitting