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:
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).
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()
fonte
Respostas:
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 é:
e no caso glm () é:
fonte
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.
fonte
R
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.
fonte
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.A escolha é baseada em sua hipótese em sua variável.
a distribuição gama é baseada em
A transformação do log baseia-se na hipótese de que,
Nesse caminho,
Com base na regra de Taylor,
Nós temos
Portanto,
No entanto, a distribuição gama baseia-se na hipótese de que,
fonte