我的问题与估计 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()
合理的?
最佳答案
使用 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 inlm
'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 的功能。
关于r - 使用 lm()、nls()(和 glm()?)估计马尔萨斯增长模型中的人口增长率,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/38001345/