glm gaussien avec fonction de lien log

droopy a écrit:

Eric Wajnberg a écrit:3) Tout ceci n'explique pas que l'usage d'une variable explicative discrète (facteur) fait que l'usage d'un lien log donne exactement les mêmes résultats qu'un lien identity.

Après réflexion je pense que c'est justement la nature discrète qui fait que l'usage d'un lien log donne les mêmes résultats (enfin presque, les coefficients ne sont pas les mêmes, mais les valeurs prédites dans l'espace de Y oui).
Pour chaque modalité, il n'y a qu'une seule valeur prédite, l'estimation de l'espérance pour cette catégorie, donc on a une constante.
Au final tu te retrouves avec E(Y|xi) = a ==> log(E(Y|xi)) = log(a), ou alors E(Y|xi) = exp(A) ==> log(E(Yxi)) = A, tu as donc que a =exp(A).
D'ailleurs l'exemple sur le facteur te montre bien que la variance est constante sur l'espace de Y puisque les modèles avec et sans lien log te donnent les mêmes estimations dans l'espace de la variable mais aussi la même estimation de la dispersion.
Il n'y a que la valeur et l'interprétation des coefficients qui changent.

Oui, je pense que c'est correct. Je vais d'abord continuer à chercher comment écrire la log-vraisemblance du modèle en question, et la décliner dans le cas d'un facteur.

Encore merci,

Eric.

niaboc a écrit:1) es-tu certain de ça? ça ressemble plus à une transformation logarithmique ce que tu écris.
Un lien log, je vois ça plutôt comme : log(E(Y))=a+bx => E(Y)=exp(a+bx) => y=exp(a+bx)+epsilon

Oui, c'est probablement correct.

Eric.

Eric Wajnberg a écrit:Pour epsilon, il me semble plus logique qu'il soit dans l'exponentiel.

Justement non, mais c'est ce que tu ne veux pas voir.
En effet un lien log veut donc dire que la linéarité va se faire sur la fonction de lien.
Donc tu as un truc du genre :

Si epsilon est dans l'espace du lien :

Quelle serait la différence entre cette formulation est celle d'un modèle ou on aurait log transformée la variable réponse avant ?

Il n'y en aurait aucune ! Ca serait le même modèle. Or quand tu appliques un lien log dans ton glm tu n'obtiens pas la même chose qu'avec un modèle sur une transformation log de Y. La différence se situe sur la variance résiduelle.
Dans le cas de la transformation log de Y tu supposes que la moyenne et la variance sont indépendants dans l'espace log (additive), ce qui implique qu'elles sont multiplicatives dans l'espace de Y.
Avec le lien log tu supposes que la variance des résidus est constante dans l'espace de Y et qu'elle est additive dans l'espace de y et non dans l'espace du lien. Donc moyenne et variance sont additives dans l'espace de Y avec ceci :

La différence entre le modèle avec le lien log et le modèle avec la variable y log transformée se situe bien au niveau de la variance.
Et ça tu le vois bien quand tu estimes les deux modèles et que tu regardes les dispersions estimées. Dans le premier cas tu as une dispersion qui correspond à la variance résiduelle sur l'espace de Y en faisant

Alors que dans le deuxième cas tu as une dispersion qui correspond à :


J'espère que ceci sera plus clair maintenant.

cdlt

Oui, oui, j'avais déjà bien compris ton explication. Tout ceci est actuellement sur le coin de mon bureau en attente (j'ai d'autres trucs urgents à régler d'abord).

Je reprendrai tout ceci prochainement, et ré-écrirai les calculs de manière exhaustive (pour moi), et les exemples numériques qui vont avec sur des jeux de données.

Je reviendrai ici si nécessaire.

Encore merci à tous pour votre aide.

Eric.

Bon, j’ai donc à présent trouvé le temps de reprendre tout ceci dans le détail, et j’arrive à présent à écrire correctement la vraisemblance d’un modèle gaussien avec un lien log, ce qui était ma question du départ et qui était restée jusqu’à présent sans réponse.

Au cas où ceci intéresserait quelqu’un ici, et comme je m’y étais engagé, je détaille la solution.

niaboc et droopy, vous aviez raison. Voici un résumé des équations, et les vérifications numériques sous R.

La variable explicative est continue

On prendra dans ce cas la variable explicative x et la variable à expliquer y suivantes :

Code:

> set.seed(1)
> x=rnorm(30 ,500, 200)
> y=rnorm(30 ,500, 200)


Commençons pas une simple régression linéaire.

Le modèle est :


avec :


La vraisemblance d’un modèle gaussien d’une simple régression linéaire dans ce cas vaut :



Et la log-vraisemblance vaut donc :



Accessoirement, si on suppose que les variances sont constantes (homoscédasticité), trouver les paramètres a et b qui maximisent cette log-vraisemblance revient à trouver les paramètres a et b qui minimisent :



C’est le critère habituel des moindres carrés.

La log-vraisemblance du modèle nul (sans la variable explicative x) vaut :



avec :



On peut recalculer tout ceci facilement à la main, sous R :

Pour le modèle nul :

Code:

> mu=mean(y)
> sd.error=sqrt(var(y)*(length(y)-1)/length(y))
> -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((y-mu)*(y-mu))
[1] -194.1389


Ce qui correspond bien au calcul fait par glm() :

Code:

> logLik(glm(y~1))
'log Lik.' -194.1389 (df=2)


Pour le modèle de régression, il faut calculer l'écart-type résiduel par rapport au modèle, qui revient - dans le cas d'un modèle avec un lien identity - à la somme des écarts des résidus du modèle par rapport à zéro, avec n-p-1 ddl, où p est le nombre de régresseurs (ici 1), soit n-2 ddl :



où les sont les résidus.

Sous R, ceci donne :

Code:

> fit=glm(y~x)
> fit$coefficients
 (Intercept)            x
504.92098572   0.04188631
> a=0.04188631
> b=504.92098572
> sd.error=sqrt(sum(fit$residuals*fit$residuals)/length(y))
>  -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((y-(a*x+b))*(y-(a*x+b)))
[1] -194.1033

Ce qui correspond encore bien au calcul fait par glm() :

Code:

> logLik(fit)
'log Lik.' -194.1033 (df=3)


Reprenons le même problème mais cette fois en log-transformant la variable y avant l’ajustement, par exemple car on soupçonne une loi log-normale sous-jacente.

Dans ce cas, le modèle est :


Ce qui revient à :



La log-vraisemblance du modèle de régression vaut cette fois-ci :



Et la log-vraisemblance du modèle nul (sans la variable explicative x) vaut :



avec cette fois-ci :



Et tout le même raisonnement suit en séquence. Sur le même jeu de données, on trouve, pour le modèle nul :

Code:

> mu=mean(log(y))
> sd.error=sqrt(var(log(y))*(length(y)-1)/length(y))
> -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((log(y)-mu)*(log(y)-mu))
[1] -8.130239


Ce qui correspond au calcul de glm() :

Code:

> logLik(glm(log(y)~1))
'log Lik.' -8.130239 (df=2)


Pour le modèle de régression :

Code:

> fit=glm(log(y)~x)
> fit$coefficients
 (Intercept)            x
6.1314714896 0.0001690827
> a=0.0001690827
> b=6.1314714896
> sd.error=sqrt(sum(fit$residuals*fit$residuals)/length(y))
>  -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((log(y)-(a*x+b))*(log(y)-(a*x+b)))
[1] -7.988918


Ce qui correspond encore bien au calcul fait par glm() :

Code:

> logLik(fit)
'log Lik.' -7.988918 (df=3)


A présent, partons plutôt sur un modèle gaussien mais avec un lien log.

Dans ce cas, le modèle est :


Commençons par le modèle nul. Il s’écrit :


Ce qui implique que :


La log-vraisemblance de ce modèle nul est :



Si on remplace :


par sa valeur :


on retombe sur la log-vraisemblance d’un modèle nul gaussien avec un lien identity :



Cette vraisemblance (pour un modèle nul) est donc la même qu’il y ait un lien log ou non, et effectivement :

Code:

> logLik(glm(y~1))
'log Lik.' -194.1389 (df=2)
> logLik(glm(y~1,family=gaussian(link="log")))
'log Lik.' -194.1389 (df=2)


Avec ce raisonnement, on peut écrire la log-vraisemblance du modèle de régression. Elle vaut :



Avec R, et dans la mesure où le lien n'est plus un lien identity, estimons l'écart-type résiduel au modèle avec la fonction sigma(), ce qui donne :

Code:

> fit=glm(y~x,family=gaussian(link="log"))
> fit$coefficients
 (Intercept)            x
6.228336e+00 7.345384e-05
> a=7.345384e-05
> b=6.228336e+00
> sd.error=sigma(fit)*sqrt((length(y)-2)/length(y))
> -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((y-exp(a*x+b))*(y-exp(a*x+b)))
[1] -194.1062
> logLik(fit)
'log Lik.' -194.1062 (df=3)


En bref, une régression sur des données log-transformées est indiquée si on modélise une variable ayant une distribution log-normale. En revanche, un glm gaussien avec un lien log permet de modéliser une variable gaussienne, mais en cherchant à ajuster une relation logarithmique entre x et y.

La variable explicative est un facteur

Passons à un schéma de type ANOVA (ou test t) (i.e., avec un lien identity), et gardons la même variable à expliquer y, mais cette fois par un facteur X à deux modalités :

Code:

> X=factor(c(rep("a",15),rep("b",15)))

Dans ce cas spécifique, la modalité “a” tient lieu par défaut de baseline codée “0”, et la modalité “b” est codée 1.

Dans ce cas, la log-vraisemblance du modèle nul reste celle d’un modèle avec un lien identity :



Comme on peut le vérifier avec R :

Code:

> sd.error=sqrt(var(y)*(length(y)-1)/length(y))
> -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((y-mean(y))*(y-mean(y)))
[1] -194.1389
> logLik(glm(y~1))
'log Lik.' -194.1389 (df=2)


Et la log-vraisemblance du modèle d’ANOVA vaut :



Sous R, ceci donne :

Code:

>> fit=glm(y~X)
> fit$coefficients
(Intercept)          Xb
   518.0572     16.9955
> a=16.9955
> b=518.0572
> sd.error=sqrt(sum(fit$residuals*fit$residuals)/length(y))
> -log(sd.error*sqrt(2*pi))*length(y)-1/(2*sd.error*sd.error)*sum((y[1:15]-b)*(y[1:15]-b))-1/(2*sd.error*sd.error)*sum((y[16:30]-(a+b))*(y[16:30]-(a+b)))
[1] -194.0945
> logLik(fit)
'log Lik.' -194.0945 (df=3)


Avec un lien log, la log-vraisemblance du modèle nul reste la même :



Et la log-vraisemblance du modèle d’ANOVA vaut :



L’ajustement du modèle donne les paramètres suivants :

Code:

> glm(y~X,family=gaussian(link="log"))$coefficients
(Intercept)          Xb
 6.25008560  0.03227958


Et on constate que l’intercept dans le modèle avec le lien log (ici 6.25008560) est le log de l’intercept avec le lien identity :

Code:

> log(518.0572)
[1] 6.250086


Ceci découle naturellement du fait que, avec le lien identity, le modèle nul est :


Alors qu’avec le lien log il est :


Il y a bien une liaison log entre les intercepts des deux modèles.

De même, (a+b) avec le lien log (ici 6.25008560+0.03227958) est bien le log de (a+b) avec le lien identity :

Code:

> 6.25008560+0.03227958
[1] 6.282365
> log(518.0572+16.9955)
[1] 6.282365


Ici aussi, ceci découle naturellement du fait qu’avec le lien identity, le modèle avec le facteur est :


Alors qu’avec le lien log il est :


Il y a bien une liaison log entre les sommes des deux coefficients des deux modèles.

En remplaçant dans les équations on retrouve bien que les deux log-vraisemblances sont bien égales avec le lien log ou le lien identity, comme on peut effectivement le calculer avec R :

Code:

>  logLik(glm(y~X))
'log Lik.' -194.0945 (df=3)
>  logLik(glm(y~X,family=gaussian(link="log")))
'log Lik.' -194.0945 (df=3)


Et tout ceci se généralise facilement avec un facteur ayant plus que deux modalités.

Tout ce raisonnement provient du fait que, au départ, j'ai trouvé contre-intuitif que - dans le cadre d'une ANOVA - l'utilisation d'un lien log ou identity ne change rien. Je crois qu'on pourrait simplement dire que, s'il existe une différence entre les moyennes des modalités d'un facteur, cette différence reste la même si on les regarde dans un espace log ou pas.

J’espère que tout ceci sera utile pour d’autres.

Eric.

Carrément.
Merci pour avoir partager tout ceci, je ne doute pas que ça soit très utile pour plein de gens.