Escolhendo entre LM e GLM para uma variável de resposta transformada em log

55

Estou tentando entender a filosofia por trás do uso de um modelo linear generalizado (GLM) vs um modelo linear (LM). Criei um exemplo de conjunto de dados abaixo em que:

registro(y)=x+ε

O exemplo não possui o erro em função da magnitude de y , portanto, eu suporia que um modelo linear de y transformado em log seria o melhor. No exemplo abaixo, esse é realmente o caso (eu acho) - já que a AIC do LM nos dados transformados em log é mais baixa. O AIC da distribuição Gamma GLM com uma função de log-link possui uma soma mais baixa de quadrados (SS), mas os graus adicionais de liberdade resultam em um AIC ligeiramente mais alto. Fiquei surpreso que a distribuição gaussiana da AIC seja muito maior (mesmo que o SS seja o mais baixo dos modelos).εy

Espero obter alguns conselhos sobre quando alguém deve abordar os modelos GLM - ou seja, há algo que eu deva procurar no meu modelo LM para ajustar os resíduos para me dizer que outra distribuição é mais apropriada? Além disso, como proceder na seleção de uma família de distribuição apropriada.

Muito obrigado antecipadamente por sua ajuda.

[EDIT]: Agora ajustei as estatísticas de resumo para que o SS do modelo linear transformado por log seja comparável aos modelos GLM com a função log-link. Um gráfico das estatísticas agora é mostrado.

Exemplo

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

insira a descrição da imagem aqui

insira a descrição da imagem aqui

Marc na caixa
fonte
exp(Xbetuma^)y1 1/2×sEugmuma2
11
Outro modelo, para o qual R não oferece uma família, é uma distribuição lognormal. SAS vai se encaixar nisso, não sei por que R glm não. Alguns sugerem gamlss do pacote R para o tgat, mas isso nunca funciona de maneira compreensível para mim. Talvez você tenha melhor sorte.
pauljohn32

Respostas:

23

Bom esforço para pensar sobre esta questão. Aqui está uma resposta incompleta, mas algumas entradas para os próximos passos.

Primeiro, as pontuações da AIC - com base nas probabilidades - estão em escalas diferentes devido às diferentes distribuições e funções de link, portanto, não são comparáveis. Sua soma dos quadrados e a soma média dos quadrados foram calculadas na escala original e, portanto, estão na mesma escala, portanto, podem ser comparadas, embora se esse seja um bom critério para a seleção de modelos seja outra questão (pode ou não ser - pesquise nos arquivos validados cruzados na seleção de modelos para uma boa discussão sobre isso).

Para sua pergunta mais geral, uma boa maneira de focar no problema é considerar a diferença entre LOG.LM (seu modelo linear com a resposta como log (y)); e LOG.GAUSS.GLM, o glm com a resposta como ye uma função de link de log. No primeiro caso, o modelo que você está ajustando é:

registro(y)=Xβ+ϵ ;

e no caso glm () é:

registro(y+ϵ)=Xβ

ϵN(0 0,σ2)

Peter Ellis
fonte
3
ϵ
4
E(Y)=g-1 1(Xβ)g(E(Y))=XβE(Y)
Achei isso muito útil: christoph-scherber.de/content/PDF%20Files/…
Aditya
16

E[em(Y|x)]em([E(Y|X]) não são os mesmos. Além disso, as suposições de variação feitas pelo GLM são mais flexíveis do que no OLS e para determinadas situações de modelagem. já que a variação da contagem pode ser diferente, tendo famílias de distribuição distintas.

Sobre a família de distribuição, na minha opinião, há uma pergunta sobre a variação e sua relação com a média. Por exemplo, em uma família gaussiana, temos variação constante. Em uma família gama, temos a variação como função quadrática da média. Plote seus resíduos padronizados versus os valores ajustados e veja como eles são.

D.Castro
fonte
11
+1 para realmente relacionada com a questão de como escolher a família certa (e eu diria que há espaço para mais alguma elaboração aqui)
etov
7

Rregistro(y)=x+εx=registro(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

As métricas para este modelo (como o AIC) não serão comparáveis ​​aos seus modelos. No entanto, sabemos que esse é o modelo certo com base no processo de geração de dados e percebemos que os coeficientes estimados estão no alvo.

- Reinstate Monica
fonte
Obrigado por seu comentário. Admito que os dados de exemplo poderiam ter sido melhores, mas acredito que estejam corretos na maneira como geraram erros. No exemplo, não há interceptação e a inclinação é 1. Se você virar a linha x = log(y) - rnorm(n, mean=0, sd=1), obterá log (y) = x + rnorm (n, média = 0, sd = 1). Se o comentário de @ whuber gerou sua resposta (acredito que sim), acredito que ele não está se referindo à geração de dados, mas à formulação do modelo GLM de @peterellis.
Marc na caixa
0

A escolha é baseada em sua hipótese em sua variável.

Vumar(XtE(Xt)=constumant

a distribuição gama é baseada em

Vumar(Xt)E(Xt)=constumant

A transformação do log baseia-se na hipótese de que,

Vumar(Xt=E(Xt)σ

Nesse caminho,

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

Com base na regra de Taylor,

registro(1 1+x)x

Nós temos

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

Portanto,

Xt=E(Xt)(1 1+Xt-E(Xt)E(Xt))registroXt=registroE(Xt)+registro(1 1+Xt-E(Xt)E(Xt))=registroE(Xt)+Xt-E(Xt)E(Xt)E(registroXt)registroE(Xt)

No entanto, a distribuição gama baseia-se na hipótese de que,

YΓ(α,β)

{E(yEu)=αEuβEuVumar(yEu)=αEuβEu2Vumar(yEu)E(yEu)=βEu
Jiaxiang
fonte