Wybór pomiędzy LM i GLM dla zmiennej odpowiedzi przekształconej logarytmicznie

55

Staram się zrozumieć filozofię stojącą za używaniem Uogólnionego Modelu Liniowego (GLM) vs Modelu Liniowego (LM). Poniżej utworzyłem przykładowy zestaw danych, w którym:

log(y)=x+ε

W przykładzie nie ma błędu w funkcji wielkości y , więc założyłbym , że najlepszy byłby model liniowy transformowanego logarytmicznie y. W poniższym przykładzie tak właśnie jest (tak mi się wydaje) - ponieważ AIC LM na danych przekształconych logami jest najniższy. AIC rozkładu gamma GLM z funkcją log-link ma niższą sumę kwadratów (SS), ale dodatkowe stopnie swobody skutkują nieco wyższym AIC. Byłem zaskoczony, że rozkład AIC Gaussa jest o wiele wyższy (mimo że SS jest najniższy z modeli).εy

Mam nadzieję uzyskać poradę, kiedy należy podejść do modeli GLM - tj. Czy jest coś, czego powinienem szukać w moim modelu LM, aby dopasować resztki, aby powiedzieć mi, że inny rozkład jest bardziej odpowiedni? Jak należy postępować w wyborze odpowiedniej rodziny dystrybucyjnej.

Z góry dziękuję za pomoc.

[EDYCJA]: Teraz dostosowałem statystyki podsumowujące, aby SS modelu liniowego przekształconego logarytmicznie było porównywalne z modelami GLM z funkcją log-link. Wykres statystyk jest teraz wyświetlany.

Przykład

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()

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

Marc w pudełku
źródło
exp(Xbeta^)y1/2×sigma2
1
Innym modelem, dla którego R nie oferuje rodziny, jest rozkład logarytmiczny. SAS to pasuje, nie wiem, dlaczego R glm nie. Niektórzy sugerują pakiet R gamlss dla tgat, ale dla mnie to nigdy nie działa zrozumiale. Może będziesz miał więcej szczęścia.
pauljohn32

Odpowiedzi:

23

Dobry wysiłek do przemyślenia tego problemu. Oto niepełna odpowiedź, ale kilka wstępów do kolejnych kroków.

Po pierwsze, wyniki AIC - oparte na prawdopodobieństwach - są w różnych skalach z powodu różnych rozkładów i funkcji łączenia, więc nie są porównywalne. Twoja suma kwadratów i średnia suma kwadratów zostały obliczone na oryginalnej skali, a zatem są na tej samej skali, więc można je porównać, chociaż to, czy jest to dobre kryterium wyboru modelu, to inne pytanie (może być, czy nie - przeszukaj sprawdzone archiwa dotyczące wyboru modelu, aby uzyskać dobrą dyskusję na ten temat).

W przypadku bardziej ogólnego pytania dobrym sposobem skoncentrowania się na problemie jest rozważenie różnicy między LOG.LM (modelem liniowym z odpowiedzią jako log (y)); i LOG.GAUSS.GLM, glm z odpowiedzią jako y i funkcją linku dziennika. W pierwszym przypadku dopasowany model to:

log(y)=Xβ+ϵ

aw przypadku glm () jest to:

log(y+ϵ)=Xβ

ϵN(0,σ2)

Peter Ellis
źródło
3
ϵ
4
E(Y)=g1(Xβ)g(E(Y))=XβE(Y)
Znalazłem to bardzo pomocne: christoph-scherber.de/content/PDF%20Files/…
Aditya
16

E[ln(Y|x)]ln([E(Y|X])

O rodzinie dystrybucyjnej moim zdaniem chodzi o wariancję i jej związek ze średnią. Na przykład w rodzinie gaussowskiej mamy stałą wariancję. W rodzinie gamma mamy wariancję jako kwadratową funkcję średniej. Wykreśl standaryzowane wartości resztkowe względem dopasowanych wartości i zobacz, jak są.

D.Castro
źródło
1
+1 za faktyczne odniesienie się do pytania o to, jak wybrać odpowiednią rodzinę (i powiedziałbym, że jest tu miejsce na trochę więcej rozwinięć)
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

Dane dla tego modelu (podobnie jak AIC) nie będą porównywalne z twoimi modelami. Wiemy jednak, że jest to właściwy model oparty na procesie generowania danych, i zauważamy, że szacowane współczynniki są właściwe dla celu.

gung - Przywróć Monikę
źródło
Dzięki za komentarz. Przyznaję, że przykładowe dane mogły być lepsze, ale uważam, że są poprawne pod względem generowania błędów. W tym przykładzie nie ma przecięcia, a nachylenie wynosi 1. Jeśli zawrócisz wokół linii x = log(y) - rnorm(n, mean=0, sd=1), otrzymasz log (y) = x + rnorm (n, średnia = 0, sd = 1). Jeśli komentarz @ Whubera zrodził twoją odpowiedź (myślę, że tak), to sądzę, że nie odnosi się on do generowania danych, ale raczej do sformułowania modelu GLM przez @peterellis.
Marc w pudełku
0

Wybór opiera się na twojej hipotezie na twojej zmiennej.

Var(XtE(Xt)=constant

rozkład gamma jest oparty na

Var(Xt)E(Xt)=constant

Transformacja logów opiera się na hipotezie, że:

Var(Xt=E(Xt)σ

W ten sposób,

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

W oparciu o zasadę Taylora,

log(1+x)x

Dostajemy

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

A zatem,

Xt=E(Xt)(1+XtE(Xt)E(Xt))logXt=logE(Xt)+log(1+XtE(Xt)E(Xt))=logE(Xt)+XtE(Xt)E(Xt)E(logXt)logE(Xt)

Jednak rozkład gamma opiera się na hipotezie, że:

YΓ(α,β)

{E(yi)=αiβiVar(yi)=αiβi2Var(yi)E(yi)=βi
Jiaxiang
źródło