Choix entre LM et GLM pour une variable de réponse transformée par un journal

55

J'essaie de comprendre la philosophie qui sous-tend l'utilisation d'un modèle linéaire généralisé (GLM) par rapport à un modèle linéaire (LM). J'ai créé un exemple de jeu de données ci-dessous où:

bûche(y)=X+ε

L'exemple n'a pas l'erreur en fonction de la magnitude de , donc je supposerais qu'un modèle linéaire de y transformé en log serait le meilleur. Dans l'exemple ci-dessous, c'est effectivement le cas (je pense), car l'AIC du LM sur les données transformées en journal est la plus basse. L'AIC du modèle GLM à distribution gamma avec fonction de liaison de journal a une somme de carrés (SS) inférieure, mais les degrés de liberté supplémentaires entraînent un AIC légèrement supérieur. J'ai été surpris de constater que la distribution gaussienne AIC est tellement plus élevée (même si le SS est le plus bas des modèles).εy

J'espère obtenir des conseils sur le moment opportun pour aborder les modèles GLM - c.-à-d. Y a-t-il quelque chose que je devrais rechercher dans les valeurs résiduelles de mon modèle LM pour me dire qu'une autre distribution est plus appropriée? En outre, comment choisir une famille de distribution appropriée?

Merci d'avance pour votre aide.

[EDIT]: J'ai maintenant ajusté les statistiques récapitulatives afin que le SS du modèle linéaire transformé en log soit comparable aux modèles GLM avec la fonction log-link. Un graphique des statistiques est maintenant affiché.

Exemple

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

entrez la description de l'image ici

entrez la description de l'image ici

Marc dans la boîte
la source
La formule de la valeur prédite dans log.lm est incorrecte. donne la médiane de . Pour obtenir la valeur attendue, ajoutez dans l'exposanty une / 2 × s i g m a 2eXp(Xbetune^)y1/2×sjegmune2
pauljohn32
1
Un autre modèle, pour lequel R n'offre pas de famille, est une distribution lognormale. SAS conviendra cela, je ne sais pas pourquoi R glm ne le fait pas. Certains suggèrent que le paquetage R soit pour tgat, mais cela ne fonctionne jamais de façon compréhensible pour moi. Peut-être aurez-vous plus de chance.
pauljohn32

Réponses:

23

Bon effort pour réfléchir à cette question. Voici une réponse incomplète, mais quelques points de départ pour les prochaines étapes.

Premièrement, les scores AIC - basés sur les probabilités - sont sur des échelles différentes en raison des distributions différentes et des fonctions de liaison différentes, et ne sont donc pas comparables. Votre somme de carrés et votre somme moyenne de carrés ont été calculées sur l’échelle initiale et sont donc sur la même échelle. Vous pouvez donc les comparer, bien qu’il s’agisse d’un bon critère pour la sélection du modèle ou non. - recherchez dans les archives croisées sur la sélection du modèle une bonne discussion à ce sujet).

Pour votre question plus générale, un bon moyen de se concentrer sur le problème consiste à considérer la différence entre LOG.LM (votre modèle linéaire avec la réponse sous la forme log (y)); et LOG.GAUSS.GLM, le gsm avec la réponse sous la forme y et une fonction de liaison de journal. Dans le premier cas, le modèle que vous montez est le suivant:

bûche(y)=Xβ+ε ;

et dans le cas glm () c'est:

bûche(y+ε)=Xβ

et dans les deux cas, est distribué .N ( 0 , σ 2 )εN(0,σ2)

Peter Ellis
la source
3
ε
4
E(Y)=g-1(Xβ)g(E(Y))=XβE(Y)
J'ai trouvé cela très utile: christoph-scherber.de/content/PDF%20Files/…
Aditya le
16

E[dans(Y|X)]dans([E(Y|X])

À mon avis, la famille de distribution est une question de variance et de sa relation avec la moyenne. Par exemple, dans une famille gaussienne, la variance est constante. Dans une famille gamma, la variance est une fonction quadratique de la moyenne. Tracez vos résidus standardisés par rapport aux valeurs ajustées et voyez comment ils sont.

D.Castro
la source
1
+1 pour avoir trait à la question de savoir comment choisir la bonne famille (et je dirais qu'il reste encore de la place pour plus de précisions ici)
jeudi
7

Rbûche(y)=X+εX=bûche(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

Les métriques pour ce modèle (comme l'AIC) ne seront pas comparables à vos modèles. Cependant, nous savons que c'est le bon modèle basé sur le processus de génération de données, et nous remarquons que les coefficients estimés sont exactement dans les objectifs.

gung - Rétablir Monica
la source
Merci pour votre commentaire. J'admets que les exemples de données auraient pu être meilleurs, mais je pense que les erreurs générées sont correctes. Dans l'exemple, il n'y a pas d'interception et la pente est 1. Si vous contournez la ligne x = log(y) - rnorm(n, mean=0, sd=1), vous obtenez log (y) = x + rnorm (n, moyenne = 0, sd = 1). Si le commentaire de @ whuber a engendré votre réponse (je crois que oui), alors je crois qu'il ne fait pas référence à la génération de données, mais plutôt à la formulation du modèle GLM par @peterellis.
Marc dans la boîte
0

Le choix est basé sur votre hypothèse sur votre variable.

Vuner(XtE(Xt)=constunent

la distribution gamma est basée sur

Vuner(Xt)E(Xt)=constunent

La transformation du journal repose sur l'hypothèse que,

Vuner(Xt=E(Xt)σ

De cette façon,

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

Basé sur la règle de Taylor,

bûche(1+X)X

On a

bûche(1+Xt-E(Xt)E(Xt))=Xt-E(Xt)E(Xt)

Ainsi,

Xt=E(Xt)(1+Xt-E(Xt)E(Xt))bûcheXt=bûcheE(Xt)+bûche(1+Xt-E(Xt)E(Xt))=bûcheE(Xt)+Xt-E(Xt)E(Xt)E(bûcheXt)bûcheE(Xt)

Cependant, la distribution gamma repose sur l'hypothèse que,

Y~Γ(α,β)

{E(yje)=αjeβjeVuner(yje)=αjeβje2Vuner(yje)E(yje)=βje
Jiaxiang
la source