使用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 

常數是已知值，而不是估計。這種常數在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「值」部分，你會得到它們是什麼更好的主意。在這些中，擬合值是$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是第一的。）但是因爲你的數據具有無噪音，你會得到使用時警告消息。

但是，就計算複雜性而言，首先採用線性模型進行變換是一個明顯的勝利。在統計建模中，變量變換爲非常普遍，所以沒有令人信服的理由拒絕使用線性模型來估計人口增長率。

作爲一個完整的圖片，我建議你嘗試所有三種方法來處理真實數據（或嘈雜的玩具數據）。估計和預測會有一些差異，但不可能非常好。

「後續跟進」

哈哈，感謝再次@Ben。對於glm()，我們也可以嘗試：

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

對於offset規格，我們可以在lm/glm使用offset參數，或者奔做的offset()功能。