[Résolu] Regression polytomique sous toutes ses coutures...

[Christophe Genolini]

Bonjour

Je cherche à faire une régression polytomique avec R mais je ne comprends pas trop les résultats qu'il me donne. Pire, il me donne parfois des choses un peu étrange, la régression polytomique étant significative la ou le drop1[Dn,.~.,test="Chisq") ne l'est pas et inversement.

Comme je risque d'avoir besoin régulièrement de ce type de régression, j'aimerai me plonger au assez profondément dans le problème et vraiment comprendre ce qui se cache derrière. Mais je n'ai pas trouvé de doc expliquant les mécanismes sous jacent. Savez-vous ou je peux en trouver une ?

En particulier, si je prends l'exemple suivant :

Code : Tout sélectionner

f<-c("O","O","O","N","O","N","O","N","O","N","O","N","O","N","O","N","O","N")
m<-c("C","C","S","MS","S","MS","MS","S","S","C","C","C","MS","S","S","MS","S","C")
m<-ordered(m,levels=c("S","C","MS"))
Dn <- data.frame(f=f,m=m)
table(Dn)

   m
f   S C MS
  N 2 3  3
  O 5 3  2

summary(polr(m~f,data=Dn))

Re-fitting to get Hessian

Call:
polr(formula = m ~ f, data = Dn)

Coefficients:
        Value Std. Error   t value
fO -0.9986798  0.8999012 -1.109766

Intercepts:
     Value   Std. Error t value
S|C  -1.0283  0.7229    -1.4223
C|MS  0.4634  0.6831     0.6783

Residual Deviance: 37.95013
AIC: 43.95013


je peux obtenir mon OR en prenant exp(-0.9986), mais comment puis-je retrouver "à la main" les coefficients de logit(P>=C)=a0 + a1 X ?

Merci de votre aide

Christophe

[Pierre Bady]

Bonjour,

La familles des modèles polytomiques est relativement vaste (ex. nominal, ordinal, + les sous-modèles). Dans R, il y a plusieurs fonctions qui permettent de construire ce type de modèle : par exemple, multinom (dans MASS) pour les variables polytomiques nominales ; polr (dans MASS) et lrm (dans Design) pour les variables polytomiques ordinales. Sinon, il existe aussi la library VGAM qui permet aussi d’analyser ces deux types de variables. Cette librairie proposée par Thomas Yee est très riche et vaut vraiment le détour :http://www.stat.auckland.ac.nz/~yee/VGAM/.

Pour retrouver les OR dans la fonction polr, c’est pas trop compliqué une fois que tu as ton modèle. Il faut utiliser les coefficients et pas oublier ce que l’on modélise : p(s), p(s)+p(c) et on en déduit p(Ms) avec la relation p(s)+p(c)+p(Ms)=1.

Code : Tout sélectionner

# suite de ton exemple

invlogit <- function(x) 1/(1+exp(-x))

p(s|f=O) = invlogit(-1.0282363+0.9986378)
p(s|f=N) = invlogit(-1.0282363)
p(s + c| f=O) = invlogit(0.4633625+0.9986378)
p(s+c| f= N) = invlogit(0.4633625)
p(c|f=O) = invlogit(0.4633625+0.9986378)-invlogit(-1.0282363+0.9986378)
p(c|f=N)= invlogit(0.4633625)-invlogit(-1.0282363)
p(Ms|f=O) = 1-invlogit(0.4633625+0.9986378)
p(Ms|f=N) = 1-invlogit(0.4633625)

# avec les fonctions lrm et predict, tu peux retrouver tous tes petits :
lrm1 <- lrm(m~f,data=Dn)
lrm1
predict(lrm1,type="fitted")
predict(lrm1,type="fitted.in")

Une fois que tu as les probabilités, tu peux calculer tes odds-ratio en utilisant la formule classique : OR= [p1/[1-p1]] / [p0/[1-p0]]. Il suffit de choisir ta référence (p0).

Par contre avec plusieurs co-variables, ça devient vite compliquer. L’estimation des intervalles de confiance devient aussi très très vite merdique. Et puis, si tu utilises des variables continues … ben … manipuler les OR avec des variables continues c’est aussi très très très délicat … il faut bien choisir la bonne échelle.

Personnellement, j’aime pas manipuler les odds-ratio et j’ai du mal à raisonner avec eux. En plus, ils sont trop souvent utiliser comme des risques relatifs. … désolé … c’est promis j’arrête de faire mon grincheux ;)


Pour la documentation, je te recommande le document pdf de Laura Thomson :

Laura Thomson’s S-PLUS (and R) Manual to Accompany Agresti’s Categorical Data Analysis. http://math.cl.uh.edu/thompsonla/Splusdiscrete2.pdf

Il existe aussi pas mal de bouquins sur le sujet (en français ou en anglais) :
- Gourieroux, C. (1989) Econométrie des variables qualitatives, Economica.
- Les bouquins d’Agresti
- etc …


en espérant t'avoir aidé un peu :)


Pierre

[Christophe Genolini]

Wahou ! Impressionnant ! Merci beaucoup.

Pour les autres fonctions, je les avais vaguement essayé suite à la lecture d'un post de ce forum ou quelqu'un en donnait déjà la liste (c'est d'ailleur ca qui m'avait fait découvrir le forum). Mais j'avais été découragé par la présentation des résultats. En effet, je travaille sur une base de 200 variables, je veux toutes les tester en bivarié et donc automatiser la récupération du petit et des OR. polr me paraissait la plus simple pour ca.

Pour ce qui qui est de retrouver les OR a la main, j'ai capté, merci. Il y a juste une bricole qui me chagrine : je pensais que la formule de base était logit(Ps)=log(Ps=1/Ps!=0)=a+bx
Or tu utilises -1.0282363+0.9986378, c'est a dire a-bx. J'esperais (sans trop y croire) que c'était une faute typographique, mais j'ai refait les calculs, c'est bien ca... D'ou sort ce changement de signe ?

En tout cas, merci beaucoup pour ton aide.

Christophe

[Pierre Bady]

bonjour,

Pour ce qui qui est de retrouver les OR a la main, j'ai capté, merci. Il y a juste une bricole qui me chagrine : je pensais que la formule de base était logit(Ps)=log(Ps=1/Ps!=0)=a+bx
Or tu utilises -1.0282363+0.9986378, c'est a dire a-bx. J'esperais (sans trop y croire) que c'était une faute typographique, mais j'ai refait les calculs, c'est bien ca... D'ou sort ce changement de signe ?

Oups … boulette ?
Dans polr, on contraint les variables explicatives à avoir un effet opposé (Agresti pp 278-279 et/ou le document pdf de Laura Thomson page 125). C'est aussi indiquer dans l'aide de polr :

Note that it is quite common for other
software to use the opposite sign for eta.

Désolé, de ne pas l’avoir préciser … :(

Si on décortique la fonction predict.polr, on peut mieux comprendre le fonctionnement du calcul des proba à partir des coefficients :

Code : Tout sélectionner

# c'est juste pour préparer les données
        newdata <- as.data.frame(Dn)
        Terms <- delete.response(plr1$terms)
        m <- model.frame(Terms, newdata, na.action = function(x) x)
        X <- model.matrix(Terms, m, contrasts = object$contrasts)
        xint <- match("(Intercept)", dimnames(X)[[2]], nomatch=0)
        if(xint > 0) X <- X[, -xint, drop=F]
        n <- nrow(X)
        q <- length(plr1$zeta)
        n
        q
        X
# coefficient de la co-variable
        plr1$coef
# valeurs des intercepts
        plr1$zeta
       
# on commence les calculs
        eta <- drop(X %*% plr1$coef)
        eta
# on voit apparaître le signe " - " :
# (le plogis correspond à invlogit du précédent post)
        cumpr <- matrix(plogis(matrix(plr1$zeta, n, q, byrow=T) - eta), , q)
        cumpr
        t(apply(cumpr, 1, function(x) diff(c(0, x, 1))))

# ouf  … on retombe sur nos pattes :)

sinon, voici l’exemple en version vglm :

Code : Tout sélectionner

require(VGAM)
vglm1 <- vglm(m~f,data=Dn,family=cumulative(parallel=T,reverse=F))
summary(vglm1)

Et pour s’amuser un peu (drôle de jeu quand même), on compare les sorties de lrm, vglm, et polr

>

Code : Tout sélectionner

 vglm1
Call:
vglm(formula = m ~ f, family = cumulative(parallel = T, reverse = F),
    data = Dn)

Coefficients:
(Intercept):1 (Intercept):2            fO
      -1.0282        0.4634        0.9986

Degrees of Freedom: 36 Total; 33 Residual
Residual Deviance: 37.95
Log-likelihood: -18.98

> plr1
Call:
polr(formula = m ~ f, data = Dn)

Coefficients:
     fO
-0.9987

Intercepts:
    S|C    C|MS
-1.0283  0.4634

Residual Deviance: 37.95
AIC: 43.95
>

Pour lrm, on a encore une autre façon de construire le modèle.
On ne modélise plus {P(Y<=1), P(Y<=2), ..., P(Y<=J-1)},
mais {P(Y>=2), P(Y>=3), ..., P(Y>=J)}.
(avec J le nombre de catégories de la variable polytomique)

Code : Tout sélectionner

> lrm1

Logistic Regression Model

lrm(formula = m ~ f, data = Dn)


Frequencies of Responses
 S  C MS
 7  6  5

       Obs  Max Deriv Model L.R.       d.f.          P          C        Dxy      Gamma      Tau-a
        18      7e-07       1.27          1     0.2607      0.607      0.215      0.418       0.15
        R2      Brier
     0.077      0.222

      Coef    S.E.   Wald Z P     
y>=C   1.0282 0.7229  1.42  0.1549
y>=MS -0.4634 0.6831 -0.68  0.4976
f=O   -0.9986 0.8999 -1.11  0.2671

> vglm(m~f,data=Dn,family=cumulative(parallel=T,reverse=T))
Call:
vglm(formula = m ~ f, family = cumulative(parallel = T, reverse = T),
    data = Dn)

Coefficients:
(Intercept):1 (Intercept):2            fO
       1.0282       -0.4634       -0.9986

Degrees of Freedom: 36 Total; 33 Residual
Residual Deviance: 37.95
Log-likelihood: -18.98
>

polr me paraissait la plus simple pour ca.

Pas sûr, …finalement ... :)


en espérant avoir éclairci mon premier post :)
@+


Pierre

[Christophe Genolini]

Pouce... Ca va trop vite pour moi... Je suis encore sur ta premiere reponse, je décortiquerais la suivante plus tard...

Pour les OR, juste pour vérification :

Code : Tout sélectionner

   m
f   S C MS
  N 2 3  3
  O 5 3  2


Grace a une regression logistique, on obtient p(s|F=N), p(s|F=O) et tout et tout. Ensuite, on doit choisir une catégorie de reférence, au hasard les personnes S (S pour sain). Puis on calcule un OR pour le C et un OR pour les MS, c'est ca ?

La formule de l'OR pour C est donc

Code : Tout sélectionner

(p(c|F=O) / p(c|F=N))  /  (p(s|F=O) / p(s|F=N))
( (2/10)  /  (3/8)  )  /  ( (5/10)  /  (2/8)  )  = 0.4


Cela sous entend qu'avoir le facteur protege, ca tombe bien, ca confirme l'intuition...

La formule de l'OR pour MS est :

Code : Tout sélectionner

(p(c|F=O) / p(c|F=N))  /  (p(s|F=O) / p(s|F=N))
( (3/10)  /  (3/8)  )  /  ( (5/10)  /  (2/8)  )  = 0.26


La regression polytomique fait ensuite l'hypothèse que l'OR est constant entre les colonnes et prend donc un OR "moyenné". Donc si on fait la moyenne entre les deux OR (0.4 et 0.26), on arrive a 0.33 c'est à dire pas tres loin du exp(-0.9986)=0.36 trouvé par la régression polytomique.

C'est ca, j'ai bon ?



Au final, je me rends compte que mes regressions polytomiques fonctionnaient... Mon problème devait venir du drop1...

Une fois la regression effectuée, pour connaitre l'effet global, j'utilise un drop1 :

Code : Tout sélectionner

polr(m~f,data=Dn)->resPol
drop1(resPol,.~.,test="Chisq")

Single term deletions

Model:
m ~ f
       Df    AIC    LRT Pr(Chi)
<none>    43.950               
f       1 43.215  1.265  0.2607


Si on compare avec un chi2 tout bete :

Code : Tout sélectionner

chisq.test(table(Dn))

        Pearson's Chi-squared test

data:  table(Dn)
X-squared = 1.2793, df = 2, p-value = 0.5275

Warning message:
l'approximation du Chi-2 est peut-être incorrecte in: chisq.test(table(Dn))


Dans un cas, p est a 0.5275, dans l'autre a 0.2607. Problème !

Je sais, je n'ai pas le droit de faire un chi2 sur des effectifs aussi petits. Mais sur des effectifs plus grands, je tombe sur le meme genre de difference...

[Pierre Bady]

bonjour,

Grace a une regression logistique, on obtient p(s|F=N), p(s|F=O) et tout et tout. Ensuite, on doit choisir une catégorie de reférence, au hasard les personnes S (S pour sain). Puis on calcule un OR pour le C et un OR pour les MS, c'est ca ?

oui, ça me semble correct.

La regression polytomique fait ensuite l'hypothèse que l'OR est constant entre les colonnes et prend donc un OR "moyenné". Donc si on fait la moyenne entre les deux OR (0.4 et 0.26), on arrive a 0.33 c'est à dire pas tres loin du exp(-0.9986)=0.36 trouvé par la régression polytomique.
C'est ca, j'ai bon ?

il semble, ... c’est ce qui est utilisé dans les autres logiciels (ex. STATA, SPSS). par contre, je n'ai pas eu le temps de le démontrer ... :(
… si c’est démontrable (à vérifier ?)

Après, je pense que le reste de tes problèmes vient du choix du type de modèle.

modèle de type 1 : logit(P(Y<=i)) = a(i) + bX
modèle de type 2 : logit(P(Y<=i)) = a(i) + b(i)X

la relation entre P(Y<=i) et X est elle commune pour tout i ?
le choix est pas nécessairement évident (cf la biblio sur le sujet).
Par défaut, polr et lrm construisent des modèles de type 1.
Avec les vglm, tu peux faire les deux.

voici un exemple qui clarifiera peut-être mieux mes propos (toujours avec les données de ton exemple initial) :

Code : Tout sélectionner

> options( prompt=" ")
# pour recuperer la librairie VGAM :
# http://www.stat.auckland.ac.nz/~yee/VGAM/download.shtml
# ou
# install.packages("VGAM", repos="http://www.stat.auckland.ac.nz/~yee")
require(VGAM)

# modèle de type 1 : logit(P(Y<=i)) = a(i) + bX

 invlogit <- function(x) 1/(1+exp(-x))
 vglm1 <- vglm(m ~ f, data=Dn, family =cumulative(parallel=T))
 coef1 <- coef(vglm1)
 coef1
(Intercept):1 (Intercept):2            fO
   -1.0282363     0.4633625     0.9986378

 psfo <-  invlogit(coef1[1]+coef1[3])
 psfn <-  invlogit(coef1[1])
 pcfo <-  invlogit(coef1[2]+coef1[3])-invlogit(coef1[1]+coef1[3])
 pcfn <-  invlogit(coef1[2])-invlogit(coef1[1])
 pmfo <-  1-invlogit(coef1[2]+coef1[3])
 pmfn <-  1-invlogit(coef1[2])
(pcfo/pcfn)/(psfo/psfn)
    0.4872274
(pmfo/pmfn)/(psfo/psfn)

    0.2605526
exp(-0.9986378)
    0.3683809
# je te l’accorde aussi,  la moyenne des deux odds-ratios est douteuse :)
 ( 0.4872274 + 0.2605526 )/2
[1] 0.37389
 
 # modèle de type 2 : logit(P(Y<=i)) = a(i) + b(i)X

 vglm2 <- vglm(m ~ f, data=Dn, family =cumulative(parallel=F))
 coef2 <- coef(vglm2)
 coef2
(Intercept):1 (Intercept):2          fO:1          fO:2
   -1.0986123     0.5108256     1.0986123     0.8754687

 psfo <-  invlogit(coef2[1]+coef2[3])
 psfn <-  invlogit(coef2[1])
 pcfo <-  invlogit(coef2[2]+coef2[4])-invlogit(coef2[1]+coef2[3])
 pcfn <-  invlogit(coef2[2])-invlogit(coef2[1])
 pmfo <-  1-invlogit(coef2[2]+coef2[4])
 pmfn <-  1-invlogit(coef2[2])
# on remarque que l’on récupère nos “bons” odds-ratios
(pcfo/pcfn)/(psfo/psfn)
          0.4
(pmfo/pmfn)/(psfo/psfn)
    0.2666667
exp(-(coef2[3]+coef2[4])/2)
    0.372678
# proche de la valeur observée pour le modèle 1 …


Le test de Chi2 est plus comparable au modèle de type 2 (valeurs de b différentes en fonction des catégories). Sa construction est plus en accord avec la structure de la table de contingence. On retrouve notamment le même nombre de degré de liberté, ce qui n’est pas le cas pour modèle 1 (cf exemple en dessous).

Code : Tout sélectionner

# pour récupérer la fonction « vanova » :
# voir dans  http://pierre.bady.free.fr/biofun/R/UtilityVglm.R
# (attention les fonctions ne sont pas toutes à jour,
# la doc est en construction et il reste certainement des bugs)

vglm0 <- vglm(m ~ 1, data=Dn, family =cumulative)
 vglm1 <- vglm(m ~ f, data=Dn, family =cumulative(parallel=T))
 vglm2 <- vglm(m ~ f, data=Dn, family =cumulative(parallel=F))
 vanova(list.vglm(list(vglm0,vglm1)),test="chisq")
Analysis of Deviance Table

Terms added sequentially (first to last)
        df     devi resdf   resdevi    pvalue
model.1 34 39.21515    34 39.215148        NA
model.2 33 37.95013     1  1.265015 0.2607039
dispersion: 1
 vanova(list.vglm(list(vglm0,vglm2)),test="chisq")

Analysis of Deviance Table

Terms added sequentially (first to last)
        df     devi resdf   resdevi    pvalue
model.1 34 39.21515    34 39.215148        NA
model.2 32 37.90819     2  1.306960 0.5202323
dispersion: 1

en espérant avoir éclairci mes précédents posts :)

@+


Pierre

PS: au fait c'est pour faire quoi ces odds-ratios ?
... désolé ... je suis curieux de nature :)

[Christophe Genolini]

Je bosse sur dépression ~ nutrition | adolescence. La question est de savoir les facteurs qui prédisposent les ado à la dépression.

Dans le milieu médical, j'ai l'impression (j'y suis depuis peu donc je ne suis pas encore sur) qu'on se fout un peu de l'effet pur, on cherche plutôt les facteurs aggravant (ou protecteur). Pour faire très bref, ce qui compte c'est de savoir que les ado ayants des tonnes de frères et soeurs sont plus protégés (OR=1.5) et que les ado faisant du sport ne sont pas spécialement protégé (après ajustement OR=1). Donc que la formule soit 10+1,5x(frereSoeur) ou 25+1,5x(frereSoeur), on s'en cogne un peu. La seule chose qui compte, c'est le 1,5 qui nous dit : les enfants uniques sont plus vulnérables à la dépression, c'est chez eux qu'il faut protéger.

Pareil, OR du sport est de 1,1 donc on arrête de gonfler les mômes qui vont pas bien pour qu'ils fassent du sport, c'est pas ca qui va les aider (contrairement aux adultes) et on pousse a ne pas sauter un repas parce que c'est clairement un facteur aggravant (OR de 2)

Voili voilou pour les OR



Concernant la régression polytomique, à priori, le modèle qu'on utilise est plutôt le modèle 1. Donc polr me va (même si c'est vachement intéressant de voir qu'on peut faire autrement. D'ailleurs, idéologiquement, quand on m'a présenté le modèle 1, j'ai trouver que c'était un peu abuser de considérer b constant...).

Par contre, je ne vois pas comment tu testes l'effet global de tes Odd ratio. Bon, on va respecter les subjects, j'ouvre un nouveau post pour le drop1.

[Pierre Bady]

rebonjour,

Dans le milieu médical, j'ai l'impression (j'y suis depuis peu donc je ne suis pas encore sur) qu'on se fout un peu de l'effet pur, on cherche plutôt les facteurs aggravant (ou protecteur)…….

bienvenu au club, … je pense comprendre ton problème.
Le milieu médical est friand d’odds-ratio … que dis-je, il en est avide … :)
Mais un odds-ratio … C’est pas évident à interpréter, à visualiser, à discuter : rapport de rapport de probabilité, … rapport des cotes,…une mesure d’association.

En ce qui concerne, les modèles polytomiques, ils nous donnent des proba, des tests de non-nullité des paramètres …toute l’information est là, que demander de plus :)

bon … j’arrête de me plaindre et de critiquer le monde médical … c’est lui qui me nourrit :)

En fait, ton « drop » marche bien, … mais dans le modèle de type 1, la co-variable a le même effet sur toutes les classes de ta variable d’intérêt, c’est la construction du modèle qu’il l’impose.

Pour rappel :
modèle de type 1 : logit(P(Y<=i)) = a(i) + bX
modèle de type 2 : logit(P(Y<=i)) = a(i) + b(i)X

Dans ton chi2, cette contrainte n’existe pas. Il peut donc exister des différences d’effets entre les classes (ce qui correspond donc plutôt au modèle de type 2). C’est pour cela que les résultats de tes deux tests divergent.


@+++

Pierre

PS: on a quand même fait un beau topic, non ? ... ;)

[Christophe Genolini]

Le milieu médical est friand d’odds-ratio … que dis-je, il en est avide … :)

Il faut dire que dans les enquetes cas témoin, seul les OR sont calculables, les RR n'ont aucun sens. Comme ils font souvent du cas témoins, ils utilisent souvent les OR et en sont familiers. Du coup, meme quand ils pourraient travailler avec autre chose, ils preferent les OR.

PS: on a quand même fait un beau topic, non ? ... ;)

Clair. D'ailleurs, je change le topic et j'y ajoute [résolu]