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

Question

我正在沿着高程梯度拟合二次模型到蜜蜂的多样性。我假设在梯度的某处有一个最大值，因此我的模型应该有一个负的“a”系数。这适用于3属，但是对于第4属（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

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任务视图，将会很有帮助。

Answer 2

我们正在处理限制估计，可以用例如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"

这暗示了线性模型。实际上，由于您的模型非常简单，因此很可能确实是理论上的最优模型，您可能需要重新考虑您的方法。例如，您是否可以将某些观察结果视为异常值并相应地改变估计值？