Elegir entre LM y GLM para una variable de respuesta transformada logarítmicamente

55

Estoy tratando de entender la filosofía detrás del uso de un modelo lineal generalizado (GLM) frente a un modelo lineal (LM). He creado un conjunto de datos de ejemplo a continuación donde:

\log (y) = x + ε

$\log(y) = x + \varepsilon$

El ejemplo no tiene el error como una función de la magnitud de , por lo que supongo que un modelo lineal de log transformado sería el mejor. En el siguiente ejemplo, este es el caso (creo), ya que el AIC del LM en los datos transformados logarítmicamente es el más bajo. El AIC del GLM de distribución Gamma con una función de enlace logarítmico tiene una suma de cuadrados (SS) más baja, pero los grados de libertad adicionales dan como resultado un AIC ligeramente más alto. Me sorprendió que la distribución gaussiana AIC sea mucho más alta (aunque el SS es el más bajo de los modelos). $\varepsilon$ $y$

Espero obtener algunos consejos sobre cuándo uno debe acercarse a los modelos GLM, es decir, ¿hay algo que deba buscar en los residuos de ajuste de mi modelo LM para decirme que otra distribución es más apropiada? Además, ¿cómo se debe proceder al seleccionar una familia de distribución adecuada?

Muchas gracias de antemano por su ayuda.

[EDITAR]: ahora he ajustado las estadísticas de resumen para que el SS del modelo lineal transformado logarítmica sea comparable a los modelos GLM con la función de enlace logarítmico. Ahora se muestra un gráfico de las estadísticas.

Ejemplo

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

ingrese la descripción de la imagen aquí

r generalized-linear-model linear-model gamma-distribution link-function Marc en la caja
fuente

La fórmula para el valor predicho en log.lm es incorrecta. da la mediana de . Para obtener el valor esperado, agregue en exponente

e x p (X \hat{b e t a})

$exp(X \hat{beta})$

y

$y$

1 / 2 \times s i g m a^{2}

$1/2 \times sigma^2$

pauljohn32

1

Otro modelo, para el cual R no ofrece una familia, es una distribución lognormal. SAS encajará en eso, no sé por qué R glm no. Algunos sugieren que el paquete R gamlss para tgat, pero nunca funciona de manera comprensible para mí. Tal vez tu tendrás mejor suerte.

pauljohn32

23

Buen esfuerzo para pensar en este tema. Aquí hay una respuesta incompleta, pero algunos entrantes para los próximos pasos.

Primero, los puntajes de AIC, basados en probabilidades, están en diferentes escalas debido a las diferentes distribuciones y funciones de enlace, por lo que no son comparables. Su suma de cuadrados y la suma media de cuadrados se han calculado en la escala original y, por lo tanto, están en la misma escala, por lo que se pueden comparar, aunque si este es un buen criterio para la selección del modelo es otra cuestión (puede ser o no) - busque en los archivos con validación cruzada sobre la selección de modelos para una buena discusión sobre esto).

Para su pregunta más general, una buena forma de enfocarse en el problema es considerar la diferencia entre LOG.LM (su modelo lineal con la respuesta como log (y)); y LOG.GAUSS.GLM, el glm con la respuesta como y y una función de enlace de registro. En el primer caso, el modelo que está ajustando es:

$\log(y)=X\beta+\epsilon$ ;

y en el caso glm () es:

$\log(y+\epsilon)=X\beta$

y en ambos casos se distribuye . $\epsilon$ $\mathcal{N}(0,\sigma^2)$

Peter Ellis
fuente

3

La caracterización de la glm no parece correcta: en el lado izquierdo hay una variable aleatoria mientras que el lado derecho contiene solo datos y parámetros, pero no variables aleatorias.

ϵ

$\epsilon$

whuber

44

Es una forma extraña de decirlo. Sé que @whuber pero proviene de convierte en . El punto es que la función de enlace va alrededor de la

E (Y) = g^{- 1} (X β)

$E(Y) = g^{-1}(X\beta)$

g (E (Y)) = X β

$g(E(Y)) = X\beta$

E (Y)

$E(Y)$

Peter Ellis

Esto me pareció muy útil: christoph-scherber.de/content/PDF%20Files/…

Aditya

16

De una manera más general, y no son lo mismo. Además, los supuestos de varianza hechos por GLM son más flexibles que en OLS, y para ciertos La situación de modelado como recuento de varianza puede ser diferente teniendo familias de distribución distintas. $E[\ln(Y|x)]$ $\ln([E(Y|X])$

Acerca de la familia de distribución en mi opinión es una pregunta sobre la varianza y su relación con la media. Por ejemplo, en una familia gaussiana tenemos una varianza constante. En una familia gamma tenemos la varianza como una función cuadrática de la media. Trace sus residuos estandarizados frente a los valores ajustados y vea cómo están.

D.Castro
fuente

1

+1 por relacionarse realmente con la cuestión de cómo elegir la familia adecuada (y diría que hay espacio para más detalles aquí)

etov

7

Desafortunadamente, su Rcódigo no conduce a un ejemplo donde . En cambio, su ejemplo es . Los errores aquí son horizontales, no verticales; son errores en , no errores en . Intuitivamente, parece que esto no debería hacer la diferencia, pero lo hace. Es posible que desee leer mi respuesta aquí: ¿Cuál es la diferencia entre la regresión lineal en y con x y x con y? Su configuración complica la cuestión de cuál es el modelo "correcto". Estrictamente, el modelo correcto es la regresión inversa: $\log(y) = x + \varepsilon$ $x = \log(y) + \varepsilon$ $x$ $y$

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

Las métricas para este modelo (como el AIC) no serán comparables con sus modelos. Sin embargo, sabemos que este es el modelo correcto basado en el proceso de generación de datos, y notamos que los coeficientes estimados están justo en el objetivo.

gung - Restablece a Monica
fuente

Gracias por tu comentario. Admito que los datos de ejemplo podrían haber sido mejores, pero creo que es correcto en cuanto a cómo generó errores. En el ejemplo, no hay intersección y la pendiente es 1. Si gira alrededor de la línea x = log(y) - rnorm(n, mean=0, sd=1), obtiene log (y) = x + rnorm (n, mean = 0, sd = 1). Si el comentario de @ whuber generó su respuesta (creo que lo hizo), entonces creo que no se está refiriendo a la generación de datos, sino a la formulación del modelo GLM de @peterellis.

Marc en la caja el

0

La elección se basa en su hipótesis sobre su variable.

la transformación del registro se basa en

\frac{\sqrt{V a r (X_{t}}}{E (X_{t})} = c o n s t a n t

$\frac{\sqrt{\mathrm{Var}(X_t} }{\mathrm{E}(X_t)} = \mathrm{constant}$

la distribución gamma se basa en

\frac{V a r (X_{t})}{E (X_{t})} = c o n s t a n t

$\frac{\mathrm{Var}(X_t) }{\mathrm{E}(X_t)} = \mathrm{constant}$

La transformación logarítmica se basa en la hipótesis de que,

\sqrt{V a r (X_{t}} = E (X_{t}) σ

$\sqrt{\mathrm{Var}(X_t} = \mathrm{E}(X_t) \sigma$

De este modo,

\begin{aligned} X_{t} & = X_{t} \\ = E (X_{t}) \cdot \frac{X_{t}}{E (X_{t})} \\ = E (X_{t}) \cdot \frac{X_{t} - E (X_{t}) + E (X_{t})}{E (X_{t})} \\ = E (X_{t}) \cdot (1 + \frac{X_{t} - E (X_{t})}{E (X_{t})}) \end{aligned}

$\begin{alignat}{2} X_t &= X_t\\ & = \mathrm{E}(X_t) \cdot \frac{X_t}{\mathrm{E}(X_t)} \\ & = \mathrm{E}(X_t) \cdot \frac{X_t - \mathrm{E}(X_t) + \mathrm{E}(X_t)}{\mathrm{E}(X_t)} \\ & = \mathrm{E}(X_t) \cdot (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ \end{alignat}$

De acuerdo con la regla de Taylor,

\log (1 + x) \approx x

$\log(1+x) \approx x$

Obtenemos

\log (1 + \frac{X_{t} - E (X_{t})}{E (X_{t})}) = \frac{X_{t} - E (X_{t})}{E (X_{t})}

$\log(1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) = \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}$

Así,

\begin{aligned} X_{t} & = E (X_{t}) \cdot (1 + \frac{X_{t} - E (X_{t})}{E (X_{t})}) \\ \log X_{t} & = \log E (X_{t}) + \log (1 + \frac{X_{t} - E (X_{t})}{E (X_{t})}) \\ = \log E (X_{t}) + \frac{X_{t} - E (X_{t})}{E (X_{t})} \\ E (\log X_{t}) & \approx \log E (X_{t}) \end{aligned}

$\begin{alignat}{2} X_t &= \mathrm{E}(X_t) \cdot (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ \log X_t &= \log \mathrm{E}(X_t) + \log (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ &= \log \mathrm{E}(X_t) + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)} \\ \mathrm{E}(\log X_t) & \approx \log \mathrm{E}(X_t) \end{alignat}$

Sin embargo, la distribución gamma se basa en la hipótesis de que,

Y \sim Γ (α, β)

$Y \sim \Gamma(\alpha, \beta)$

{\begin{cases} E (y_{i}) = α_{i} \cdot β_{i} \\ V a r (y_{i}) = α_{i} \cdot β_{i}^{2} \end{cases} \to \frac{V a r (y_{i})}{E (y_{i})} = β_{i}

$\begin{cases} E(y_i) = \alpha_i \cdot \beta_i \\ Var(y_i) = \alpha_i \cdot \beta_i^2 \\ \end{cases} \to \frac{Var(y_i)}{E(y_i)} = \beta_i$

Jiaxiang
fuente

Elegir entre LM y GLM para una variable de respuesta transformada logarítmicamente

Ejemplo

Respuestas: