Mit lm(), nls() (und glm()?) Zur Schätzung der Bevölkerungswachstumsrate im Malthusianischen Wachstumsmodell

Question

Meine Frage bezieht sich auf die Schätzung der Bevölkerungswachstumsrate in Malthusian growth model. Als Spielzeug Beispiel einen Spielzeug-Datensatz 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") 

Ich versuche, diesen Datensatz durch exponentielles Modell zu passen:

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

und r schätzen. Wenn nichtlineare Regression unter Verwendungnls():

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

bekomme ich folgende Fehlermeldung:

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

Ich habe auch versucht lm()

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

aber erhalten

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

Wie kann ich das lösen? Wie kann ich die Daten an das Exponentialmodell anpassen?

Gibt es auch andere Ansätze, die ich für die Anpassung des Wachstumsmodells in Betracht ziehen könnte? Ist glm() sinnvoll?

Answer 1

Mit lm()

Bitte lesen Sie ?formula für die korrekte Angabe einer Formel. Jetzt gehe ich davon aus, dass du das gelesen hast.

Zuerst Ihr Modell, nach der Einnahme von log sowohl auf LHS und RHS-Transformation wird:

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

Die Konstante ist ein bekannter Wert, nicht geschätzt werden. Eine solche Konstante heißt offset in lm.

Sie sollten lm wie diese verwenden:

# "-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 

Wie Sie entdeckt haben, fit ist eine Liste der Länge 13. Siehe den „Wert“ in ?lm und Sie werden bessere Vorstellung davon, was sie sind . Unter diesen sind die angepassten Werte $fitted, so können Sie Ihr Grundstück ziehen durch:

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

Achten Sie auf meine exp(fit$fitted) verwenden, weil wir ein Modell für log(y) passen und jetzt gehen wir zurück zu ursprüngliche Skala.

Bemerkung

Wie @BenBolker sagte, eine einfachere Beschreibung ist:

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

oder

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

Aber die Antwortvariable ist nicht log(y) aber log(y/10000) jetzt, so, wenn Sie Machen Sie Grundstück, Sie benötigen:

---

Verwendung nls()

Correct Weise nls() für die Verwendung ist wie folgt:

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

Da nichtlineare Kurvenanpassung Iterationen erfordert, wird ein Startwert benötigt wird, und must über das Argument start bereitgestellt werden. Jetzt

, wenn Sie diesen Code versuchen, erhalten Sie:

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

Das Problem ist, weil Ihre Daten genau sind, ohne Lärm. Haben Sie einen Read auf ?nls:

Warning: 

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

So setzen nls() für Ihre Spielzeug Daten df funktioniert nicht.

Gehen wir zurück das angepasste Modell von lm() zu überprüfen:

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 

Residuen sind 0 im Grunde überall, und lm() tut passgenau in diesem Fall.

---

Follow-up

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.

Es gibt tatsächlich einen Unterschied in der Formel zwischen lm und nls. Vielleicht können Sie es als solche nehmen:

• lm() ‚s Formel Modellformel genannt, die Sie von ?formula lesen kann. Es ist so fundamental in R. Modell Anpassungsroutinen es verwenden, wie lm, glm, während viele Funktionen Formel Verfahren haben, wie model.matrix, aggregate, boxplot usw.
• nls() ‚s Formel ist eher wie eine Funktionsspezifikation, und wirklich nicht weit verbreitet. Viele andere Funktionen, die nicht-lineare Iterationen wie optim ausführen, akzeptieren keine Formel, sondern nehmen direkt eine Funktion an. Also, behandeln Sie einfach nls() als Sonderfall.

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.

Streng genommen reale Bevölkerungsdaten geben (auf jeden Fall mit Rauschen), nls() für die Verwendung von Kurvenanpassung, oder mit glm(, family = poisson) für eine Poisson Reaktion GLM besseren Boden hat, als ein lineares Modell paßt. Der glm() Anruf auf Ihre Daten wäre:

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

(Sie müssen möglicherweise lernen, was ein GLM ist zuerst.) Aber da Ihre Daten keinen Lärm haben, werden Sie eine Warnmeldung erhalten, wenn es zu benutzen.

In Bezug auf Rechenkomplexität ist die Verwendung eines linearen Modells durch die erste Transformation log ein klarer Gewinn. In der statistischen Modellierung sind die Variablentransformation sehr häufig, es gibt also keinen zwingenden Grund, die Verwendung des linearen Modells zur Schätzung der Populationswachstumsrate abzulehnen.

Als vollständiges Bild empfehle ich Ihnen, alle drei Ansätze für echte Daten (oder laute Spielzeugdaten) zu versuchen. Es wird einige Unterschiede in der Schätzung und Vorhersage geben, aber wahrscheinlich nicht sehr groß sein.

"Follow-Follow-up"

Haha, nochmals vielen Dank an @ Ben. Für glm(), können wir auch versuchen:

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

für offset Spezifikation, entweder wir offset Argument in lm/glm oder offset() Funktion als Ben tut verwenden können.