Auswahl zwischen LM und GLM für eine log-transformierte Antwortvariable

55

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:

Log(y)=X+ε

Das Beispiel hat nicht den Fehler als Funktion der Größe vonyεy , 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()

Bildbeschreibung hier eingeben

Bildbeschreibung hier eingeben

Marc in der Kiste
quelle
exp(Xbeta^)y1/2×sigma2
1
Ein weiteres Modell, für das R keine Familie anbietet, ist eine logarithmische Normalverteilung. SAS wird das passen, ich weiß nicht, warum R glm nicht. Einige schlagen R-Paket-Gamlss für Tgat vor, aber es funktioniert nie verständlich für mich. Vielleicht haben Sie besseres Glück.
Pauljohn32

Antworten:

23

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:

Log(y)=Xβ+ϵ

und im glm () Fall ist es:

Log(y+ϵ)=Xβ

ϵN(0,σ2)

Peter Ellis
quelle
3
ϵ
4
E(Y.)=G-1(Xβ)G(E(Y.))=XβE(Y.)
Ich fand das sehr hilfreich: christoph-scherber.de/content/PDF%20Files/…
Aditya
16

E[ln(Y.|X)]ln([E(Y.|X]) nicht identisch, und die von GLM gemachten Varianzannahmen sind flexibler als in OLS und für bestimmte Modellierungssituationen Die Varianz kann bei unterschiedlichen Verteilungsfamilien unterschiedlich sein.

Ü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.

D. Castro
quelle
1
+1 für den tatsächlichen Bezug auf die Frage, wie die richtige Familie zu wählen (und ich würde sagen, es gibt Raum für weitere Ausarbeitung hier)
etov
7

RLog(y)=X+εX=Log(y)+εXy

ly = log(y)
REVERSE.REGRESSION = lm(x~ly)
summary(REVERSE.REGRESSION)
# Call:
# lm(formula = x ~ ly)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -2.93996 -0.64547 -0.01351  0.63133  2.92991 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.01563    0.03113   0.502    0.616    
# ly           1.01519    0.03138  32.350   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.984 on 998 degrees of freedom
# Multiple R-squared:  0.5119,    Adjusted R-squared:  0.5114 
# F-statistic:  1047 on 1 and 998 DF,  p-value: < 2.2e-16

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.

gung - Wiedereinsetzung von Monica
quelle
Vielen Dank für Ihren Kommentar. Ich gebe zu, die Beispieldaten hätten besser sein können, aber ich glaube, dass sie in Bezug auf die Fehlererzeugung korrekt sind. In diesem Beispiel gibt es keinen Schnittpunkt und die Steigung ist 1. Wenn Sie sich um die Linie drehen 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.
Marc in der Box
0

Die Auswahl basiert auf Ihrer Hypothese zu Ihrer Variablen.

Veinr(XtE(Xt)=cOnsteinnt

Die Gamma-Verteilung basiert auf

Veinr(Xt)E(Xt)=cOnsteinnt

Die log-Transformation beruht auf der Hypothese, dass

Veinr(Xt=E(Xt)σ

Auf diese Weise,

Xt=Xt=E(Xt)XtE(Xt)=E(Xt)Xt-E(Xt)+E(Xt)E(Xt)=E(Xt)(1+Xt-E(Xt)E(Xt))

Basierend auf der Taylor-Regel

Log(1+X)X

Wir bekommen

Log(1+Xt-E(Xt)E(Xt))=Xt-E(Xt)E(Xt)

Somit,

Xt=E(Xt)(1+Xt-E(Xt)E(Xt))LogXt=LogE(Xt)+Log(1+Xt-E(Xt)E(Xt))=LogE(Xt)+Xt-E(Xt)E(Xt)E(LogXt)LogE(Xt)

Die Gamma-Verteilung beruht jedoch auf der Hypothese, dass

Y.Γ(α,β)

{E(yich)=αichβichVeinr(yich)=αichβich2Veinr(yich)E(yich)=βich
Jiaxiang
quelle