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

Question

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

隔离这个属，我们可以清楚地看到为什么它是“不正确的”：

有一个二次模型和一个线性模型。考虑到数据点，二次方是有意义的，但在生物学上没有多大意义。我想强制命令生成负“a”（因此给出的“最佳”高度可能远低于第一张图中给出的高度，即 1193 m），我该怎么做？ 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 列出了几个解决最小二乘问题的包，允许对系数进行（线性）约束，例如 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更可靠，并尽量避免“零梯度”。
如果用户更频繁地利用许多 CRAN 任务视图，将会很有帮助。