r - 使用 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正确指定公式。现在我将继续假设您已经阅读了该内容。

首先是你的模型，拍下后log LHS 和 RHS 上的变换，变为:

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

常数是已知值，不可估计。这样的常数被称为 offset在 lm .

您应该使用 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 的“值”部分你会更好地了解它们是什么。其中，拟合值为 $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)现在，所以当你制作情节时，你需要:

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转型是一个明显的胜利。在统计建模中，变量变换是 很常见 ，因此没有令人信服的理由拒绝使用线性模型来估计人口增长率。

作为一个完整的图片，我建议您针对真实数据(或嘈杂的玩具数据)尝试所有三种方法。估计和预测会有一些差异，但不太可能很大。

《跟进》

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

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

对于 offset规范，我们可以使用 offset参数在 lm/ glm ，或 offset() Ben 的功能。