Используя lm (), nls () (и 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. См. раздел "значение" в разделе И вы получите лучшее представление о том, что они собой представляют. Среди них подходящие значения $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() в этом случае точно соответствует.

---

Последующие действия

Еще одна вещь, которую я не смог выяснить, - это почему параметр r не используется в спецификации формулы lm.

На самом деле есть некоторая разница в Формуле между lm и nls. Возможно, ты сможешь. примите это как таковое:

• Формула lm() называется формулой модели, которую можно прочитать из ?formula. Она настолько фундаментальна в R. процедуры подгонки модели используют ее, как lm, glm, в то время как многие функции имеют метод формулы, как model.matrix, aggregate, boxplot, и т.д.
• Формула nls() больше похожа на спецификацию функции и действительно не широко используется. Многие другие функции, выполняющие нелинейные итерации, такие как optim, не принимают формулу, но принимают функцию непосредственно. Так что, просто рассматривайте nls() как частный случай.

Так есть ли смысл делать это с помощью линейной модели? Просто то, что я пытаюсь смоделировать здесь, - это использование мальтузианской модели роста.

Строго говоря, предоставление реальных данных о населении (конечно, с шумом), использование nls() для подгонки кривой или использование glm(, family = poisson) для реакции Пуассона GLM имеет лучшую основу, чем подгонка линейной модели. Вызов glm() для ваших данных будет следующим:

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

(возможно, вам нужно узнать, что ГЛМ-это первое.) Но так как ваши данные не имеют шума, вы получите предупреждающее сообщение при его использовании.

Однако с точки зрения вычислительной сложности использование линейной модели путем первого преобразования logявляется очевидным выигрышем. В статистическом моделировании переменные преобразования очень распространены, поэтому нет никаких веских причин отказываться от использования линейной модели для оценки темпов роста населения.

В качестве полной картины я рекомендую вам попробовать все три подхода по-настоящему данные (или шумные игрушечные данные). Будет некоторая разница в оценках и прогнозах, но вряд ли она будет очень большой.

"Follow-follow-up"

Ха-ха, еще раз спасибо @Ben. Для glm() Мы также можем попробовать:

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

Для спецификации offset мы можем либо использовать аргумент offset в lm/glm, или функцию offset(), как это делает Бен.