Ich versuche die Philosophie zu verstehen, die hinter der Verwendung eines generalisierten linearen Modells (GLM) gegenüber einem linearen Modell (LM) steckt. Ich habe unten einen Beispieldatensatz erstellt:
Das Beispiel hat nicht den Fehler als Funktion der Größe vony , daher würde ich annehmen, dass ein lineares Modell des logarithmisch transformierten y das beste wäre. Im folgenden Beispiel ist dies tatsächlich der Fall (glaube ich), da der AIC des LM auf den log-transformierten Daten am niedrigsten ist. Der AIC der Gamma-Verteilung GLM mit einer Log-Link-Funktion hat eine geringere Quadratsumme (SS), aber die zusätzlichen Freiheitsgrade führen zu einem etwas höheren AIC. Ich war überrascht, dass der AIC der Gaußschen Verteilung so viel höher ist (obwohl der SS das niedrigste Modell ist).
Ich hoffe, einige Ratschläge zu erhalten, wann man sich GLM-Modellen nähern sollte - dh gibt es etwas, nach dem ich in meinen Residuen für LM-Modelle Ausschau halten sollte, um mir mitzuteilen, dass eine andere Verteilung angemessener ist? Wie ist bei der Auswahl einer geeigneten Distributionsfamilie vorzugehen?
Vielen Dank im Voraus für Ihre Hilfe.
[EDIT]: Ich habe jetzt die Auswertungsstatistik so angepasst, dass die SS des log-transformierten linearen Modells mit den GLM-Modellen mit der Log-Link-Funktion vergleichbar ist. Ein Diagramm der Statistik wird jetzt angezeigt.
Beispiel
set.seed(1111)
n <- 1000
y <- rnorm(n, mean=0, sd=1)
y <- exp(y)
hist(y, n=20)
hist(log(y), n=20)
x <- log(y) - rnorm(n, mean=0, sd=1)
hist(x, n=20)
df <- data.frame(y=y, x=x)
df2 <- data.frame(x=seq(from=min(df$x), to=max(df$x),,100))
#models
mod.name <- "LM"
assign(mod.name, lm(y ~ x, df))
summary(get(mod.name))
plot(y ~ x, df)
lines(predict(get(mod.name), newdata=df2) ~ df2$x, col=2)
mod.name <- "LOG.LM"
assign(mod.name, lm(log(y) ~ x, df))
summary(get(mod.name))
plot(y ~ x, df)
lines(exp(predict(get(mod.name), newdata=df2)) ~ df2$x, col=2)
mod.name <- "LOG.GAUSS.GLM"
assign(mod.name, glm(y ~ x, df, family=gaussian(link="log")))
summary(get(mod.name))
plot(y ~ x, df)
lines(predict(get(mod.name), newdata=df2, type="response") ~ df2$x, col=2)
mod.name <- "LOG.GAMMA.GLM"
assign(mod.name, glm(y ~ x, df, family=Gamma(link="log")))
summary(get(mod.name))
plot(y ~ x, df)
lines(predict(get(mod.name), newdata=df2, type="response") ~ df2$x, col=2)
#Results
model.names <- list("LM", "LOG.LM", "LOG.GAUSS.GLM", "LOG.GAMMA.GLM")
plot(y ~ x, df, log="y", pch=".", cex=3, col=8)
lines(predict(LM, newdata=df2) ~ df2$x, col=1, lwd=2)
lines(exp(predict(LOG.LM, newdata=df2)) ~ df2$x, col=2, lwd=2)
lines(predict(LOG.GAUSS.GLM, newdata=df2, type="response") ~ df2$x, col=3, lwd=2)
lines(predict(LOG.GAMMA.GLM, newdata=df2, type="response") ~ df2$x, col=4, lwd=2)
legend("topleft", legend=model.names, col=1:4, lwd=2, bty="n")
res.AIC <- as.matrix(
data.frame(
LM=AIC(LM),
LOG.LM=AIC(LOG.LM),
LOG.GAUSS.GLM=AIC(LOG.GAUSS.GLM),
LOG.GAMMA.GLM=AIC(LOG.GAMMA.GLM)
)
)
res.SS <- as.matrix(
data.frame(
LM=sum((predict(LM)-y)^2),
LOG.LM=sum((exp(predict(LOG.LM))-y)^2),
LOG.GAUSS.GLM=sum((predict(LOG.GAUSS.GLM, type="response")-y)^2),
LOG.GAMMA.GLM=sum((predict(LOG.GAMMA.GLM, type="response")-y)^2)
)
)
res.RMS <- as.matrix(
data.frame(
LM=sqrt(mean((predict(LM)-y)^2)),
LOG.LM=sqrt(mean((exp(predict(LOG.LM))-y)^2)),
LOG.GAUSS.GLM=sqrt(mean((predict(LOG.GAUSS.GLM, type="response")-y)^2)),
LOG.GAMMA.GLM=sqrt(mean((predict(LOG.GAMMA.GLM, type="response")-y)^2))
)
)
png("stats.png", height=7, width=10, units="in", res=300)
#x11(height=7, width=10)
par(mar=c(10,5,2,1), mfcol=c(1,3), cex=1, ps=12)
barplot(res.AIC, main="AIC", las=2)
barplot(res.SS, main="SS", las=2)
barplot(res.RMS, main="RMS", las=2)
dev.off()
quelle
Antworten:
Gute Mühe zum Durchdenken dieses Themas. Hier ist eine unvollständige Antwort, aber einige Vorspeisen für die nächsten Schritte.
Erstens sind die AIC-Scores - basierend auf den Wahrscheinlichkeiten - aufgrund der unterschiedlichen Verteilungen und Verknüpfungsfunktionen unterschiedlich skaliert und daher nicht vergleichbar. Die Summe der Quadrate und die mittlere Summe der Quadrate wurden auf der ursprünglichen Skala berechnet und liegen daher auf derselben Skala. Sie können also vergleichen, auch wenn dies ein gutes Kriterium für die Modellauswahl ist - Durchsuchen Sie die kreuzvalidierten Archive nach Modellauswahl, um dies zu erörtern.
Für Ihre allgemeinere Frage besteht eine gute Möglichkeit, sich auf das Problem zu konzentrieren, darin, den Unterschied zwischen LOG.LM (Ihrem linearen Modell mit der Antwort als log (y)) zu betrachten. und LOG.GAUSS.GLM, die glm mit der Antwort als y und einer Protokollverbindungsfunktion. Im ersten Fall ist das Modell, das Sie anpassen:
und im glm () Fall ist es:
quelle
Über die Verbreitungsfamilie geht es meiner Meinung nach um die Varianz und deren Beziehung zum Mittelwert. Zum Beispiel haben wir in einer Gaußschen Familie eine konstante Varianz. In einer Gammafamilie haben wir die Varianz als quadratische Funktion des Mittelwerts. Zeichnen Sie Ihre standardisierten Residuen gegen die angepassten Werte und sehen Sie, wie sie sind.
quelle
R
Metriken für dieses Modell (wie das AIC) sind nicht mit Ihren Modellen vergleichbar. Wir wissen jedoch, dass dies das richtige Modell ist, das auf dem Datenerzeugungsprozess basiert, und stellen fest, dass die geschätzten Koeffizienten genau auf dem Ziel liegen.
quelle
x = log(y) - rnorm(n, mean=0, sd=1)
, erhalten Sie log (y) = x + rnorm (n, Mittelwert = 0, sd = 1). Wenn der Kommentar von @ whuber Ihre Antwort hervorgebracht hat (ich glaube das hat er), dann bezieht er sich meiner Meinung nach nicht auf die Datengenerierung, sondern auf die GLM-Modellformulierung von @peterellis.Die Auswahl basiert auf Ihrer Hypothese zu Ihrer Variablen.
Die Gamma-Verteilung basiert auf
Die log-Transformation beruht auf der Hypothese, dass
Auf diese Weise,
Basierend auf der Taylor-Regel
Wir bekommen
Somit,
Die Gamma-Verteilung beruht jedoch auf der Hypothese, dass
quelle