使用 lm()、nls()（和 glm()？）估计马尔萨斯增长模型中的人口增长率

Question

我的问题与估计 Malthusian growth model 的人口增长率有关。作为玩具示例，考虑一个玩具数据集 df:

structure(list(x= c(0L, 24L, 48L, 72L, 96L, 120L, 144L, 168L
), y = c(10000, 18744.0760659189, 35134.0387564953, 65855.509495469, 
123440.067934292, 231377.002294256, 433694.813090781, 812920.856596808
)), .Names = c("x", "y"), row.names = c(NA, -8L), class = "data.frame")

我正在尝试通过 指数模型:

来拟合这个数据集

y = 10000 * (e^(r * x))

并估计 r。当使用 非线性 回归时 nls():

fit <- nls(y ~ (10000 * exp(r*x)), data=df)

我收到以下错误：

Error in getInitial.default(func, data, mCall = as.list(match.call(func,  : 
  no 'getInitial' method found for "function" objects

我也试过了lm()

fit <- lm(log(y) ~ (10000 * exp(r*x)), data=df)

但得到

Error in terms.formula(formula, data = data) : 
  invalid model formula in ExtractVars

我该如何解决这个问题？如何将数据拟合到我拥有的指数模型？

此外，我可以考虑使用其他方法来拟合人口增长模型吗？ glm()合理吗？

Answer 1

使用 lm()

请阅读 ?formula 以正确指定公式。现在我将继续假设你已经阅读了。

首先，您的模型在对 LHS 和 RHS 进行 log 变换后变为：

log(y) = log(10000) + r * x

该常数是已知值，不可估计。这样的常量在lm中叫做offset。

你应该这样使用 lm:

# "-1" in the formula will drop intercept
fit <- lm(log(y) ~ x - 1, data = df, offset = rep(log(10000), nrow(df)))

# Call:
#  lm(formula = log(y) ~ x - 1, data = df, offset = rep(log(10000), nrow(df)))

#  Coefficients:
#        x  
#  0.02618

如您所见，fit 是一个长度为 13 的列表。请参阅 ?lm 的 "Value" 部分，您将更好地了解它们是什么。其中，拟合值为 $fitted，因此您可以通过以下方式绘制绘图：

plot(df)
lines(df$x, exp(fit$fitted), col = 2, lwd = 2)  ## red line

注意我使用 exp(fit$fitted)，因为我们为 log(y) 拟合了一个模型，现在我们要回到原来的比例。

备注

正如@BenBolker所说，一个更简单的规范是：

fit <- lm(log(y/10000) ~ x - 1, data = df)

或

fit <- lm(log(y) - log(10000) ~ x - 1, data = df)

但是现在的响应变量不是log(y)而是log(y/10000)，所以做plot的时候需要：

lines(df$x, 10000 * exp(fit$fitted), col = 2, lwd = 2)

使用nls()

nls()的正确使用方法是这样的：

nls(y ~ 10000 * exp(r * x), data = df, start = list(r = 0.1))

因为非线性曲线拟合需要迭代，所以需要一个起始值，并且必须通过参数start.

提供

现在，如果您尝试此代码，您将得到：

Error in nls(y ~ 10000 * exp(r * x), data = df, start = list(r = 0.1)) : 
  number of iterations exceeded maximum of 50

问题是因为你的数据是准确的，没有噪音。阅读 ?nls:

Warning:

     *Do not use ‘nls’ on artificial "zero-residual" data.*

因此，将 nls() 用于您的玩具数据集 df 不起作用。

让我们返回检查来自 lm() 的拟合模型：

fit$residuals
#            1             2             3             4             5 
#-2.793991e-16 -1.145239e-16 -2.005405e-15 -5.498411e-16  3.094618e-15 
#            6             7             8 
# 1.410007e-15 -1.099682e-15 -1.007937e-15

残差基本上到处都是 0，lm() 在这种情况下完全适合。

跟进

One last thing that I haven't been able to figure out is why the parameter r is not used in lm's formula specification.

lm和nls的公式其实有些区别。或许你可以这样理解：

lm()的公式称为模型公式，可参考?formula。它在 R 中非常基础。模型拟合例程使用它，如 lm、glm，而许多函数都有公式方法，如 model.matrix、aggregate、boxplot,等等
nls()的公式更像是一个函数规范，并没有被广泛使用。许多其他进行非线性迭代的函数，如 optim 将不接受公式，而是直接采用函数。所以，就把nls()当作一个特例吧。

So would it make sense to do it using the linear model? Simply what I am trying to model here is using Malthusian growth model.

严格来说，给出真实的人口数据（当然有噪声），使用 nls() 进行曲线拟合，或使用 glm(, family = poisson) 进行泊松响应 GLM 比拟合线性模型具有更好的基础。 glm() 调用您的数据将是：

glm(y ~ x - 1, family = poisson(), data = df, offset = rep(log(10000), nrow(df)))

（您可能需要先了解什么是 GLM。）但是由于您的数据没有噪声，因此在使用它时会收到警告消息。

但是，就计算复杂度而言，首先采用 log 变换使用线性模型显然是一个胜利。在统计建模中，变量变换非常普遍，因此没有令人信服的理由拒绝使用线性模型来估计人口增长率。

总的来说，我建议您对真实数据（或嘈杂的玩具数据）尝试所有三种方法。估计和预测会有一些差异，但不会很大。

"Follow-follow-up"

哈哈，再次感谢@Ben。对于glm()，我们也可以试试：

glm(y ~ x - 1 + offset(log(10000)), family = gaussian(link="log"))

对于offset规范，我们可以在lm/glm中使用offset参数，或者像Ben那样使用offset()函数。

使用 lm()、nls()（和 glm()？）估计马尔萨斯增长模型中的人口增长率

Using lm(), nls() (and glm()?) to estimate population growth rate in Malthusian growth model

r

regression

glm

lm

nls