Régression logistique multinomiale

Bonjour à tous!

Je suis étudiant en master de biologie. Je suis en train d'écrire mon mémoire et aujourd'hui je suis face à un mur.
Cela fait un mois maintenant que j'essaye de résoudre mon problème de statistique en vain.
Voilà ou j'en suis aujourd'hui :

J'ai essayé trois programmes: R, JMP et SPSS.
J'ai abandonné JMP car je ne suis pas sur que ce soit le programme le plus approprié pour l'analyse que j'essaye de faire. Je me suis donc cantonné à R et SPSS.

Je vais en premier lieu vous décrire mes données :

J'ai une table de donnée (qui vient d'excel mais que j'ai transformé en .txt pour une importation plus facile sous R) avec (chaque point représente une variable de ma table, le titre d'une colonne):

• Date (date de l'observation)
  
• Temperature (au moment de l'observation)
  
• Behavior (comportement observé, 3 choix différents : F (Feeding), R(Resting), M(Moving), RUN(Runing))
  
• Valley (vallée ou le comportement à été observé, 2 choix possible (Fuorn et Trupchun)
  
• Age (de l'individu)
  
• Mois (de l'observation)
  
• Année (de l'osbservation)
  
• Kid (si l'individu à un enfant ou pas)
  
• Individu (numéro du tag de la base de donnée du parc national ou j'ai pris mes données, j'ai 11 individus différents)

Mon but ici est de déterminer quels facteurs influencent le comportement des individus.

Passons à présent aux analyse que j'ai faites :

Sous R

Package : VGAM
formule : mod <- vglm(Behavior ~ Temp + Valley + Age + Month + Year + Individual, family=multinomial(), data=Merge)
Summary(mod)

Résultat :

Code:

Pearson Residuals:
                      Min      1Q  Median      3Q    Max
log(mu[,1]/mu[,4]) -29.903 -0.66948  0.61220  0.82061 1.2109
log(mu[,2]/mu[,4]) -29.054 -0.27864 -0.23503 -0.20717 3.8363
log(mu[,3]/mu[,4]) -29.453 -0.54559 -0.42744  0.95586 2.7680

Code:

Coefficients:
                    Estimate Std. Error  z value
(Intercept):1    -3.2212e+02 6.8413e+01  -4.7085
(Intercept):2    -1.7247e+02 6.9163e+01  -2.4937
(Intercept):3    -4.4037e+02 6.8467e+01  -6.4318
Temp:1            4.2142e-02 4.4678e-03  9.4324
Temp:2            6.2538e-02 4.5311e-03  13.8020
Temp:3          -5.4563e-03 4.4795e-03  -1.2181
ValleyTrupchun:1 -3.8754e+00 3.6403e-01 -10.6460
ValleyTrupchun:2 -4.8341e+00 3.7123e-01 -13.0217
ValleyTrupchun:3 -5.1426e+00 3.6443e-01 -14.1113
Age:1            -1.6646e-01 1.4318e-02 -11.6256
Age:2            -1.6630e-01 1.4555e-02 -11.4262
Age:3            -1.8485e-01 1.4342e-02 -12.8886
Month:1          -9.2571e-01 3.6144e-02 -25.6117
Month:2          -9.7004e-01 3.6720e-02 -26.4174
Month:3          -8.8669e-01 3.6186e-02 -24.5036
Year:1            1.6571e-01 3.4046e-02  4.8673
Year:2            9.0428e-02 3.4419e-02  2.6273
Year:3            2.2459e-01 3.4073e-02  6.5916
Individual:1      5.4454e-03 7.1586e-04  7.6068
Individual:2      7.4979e-03 7.3060e-04  10.2626
Individual:3      8.5069e-03 7.1665e-04  11.8704

Number of linear predictors:  3

Names of linear predictors:
log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]), log(mu[,3]/mu[,4])

Dispersion Parameter for multinomial family:  1

Residual deviance: 530182.6 on 815760 degrees of freedom

Log-likelihood: -265091.3 on 815760 degrees of freedom

Number of iterations: 8

donc les soucis que je rencontre sous R sont que je ne sais pas à quoi 1, 2, 3 (ex: Kids:1, Kids:2, Kids:3) font référence ni quel est le comportement de référence. Egalement je ne sais pas comment afficher mes p-values!

J'espère que quelqu'un ici connait un peu le package VGAM pour pouvoir m'aider car R est le programme utilisé principalement par mon université.

Sous SPSS

Dans ce programme, j'ai importé mon fichier .txt et mit un label sur chacune de mes variable :



Ma table de donnée :



J'ai ensuite lancé mon analyse multinomiale :



Voici les résultats :





Donc ce qui me pose problème ici :

Je trouve une p-values très faibles pour des variables qui n'ont pas beaucoup d'effet (normalement la vallée à très peu d'effet).
Dans la table "qualité d'ajustement" je trouve des p-value très faibles alors que pour une bonne qualité la p-value doit être élevée (si j'ai bien compris).
Dans la table "Estimation des paramètres" presque tous mes paramètres ont une p-value de 0.000 également, cela me parait bizarre.

Dernièrement dans la table "classification" je trouve que ma table est a 52% correctement estimée.

En gros j'ai plein de chose à dire dessus mais il me semble qu'il y a une faute quelque part, soit dans mon data set (peut-être que mes variables sont mal agencées) soit dans l'analyse ou je n'ai pas coché les bonnes options.

La date de rédition de mon mémoire approche à grands pas et je suis complètement bloqué.
J'espère trouver de l'aide ici!

En attendant, bonne fêtes à tous!

Lapive

Personnellement j'aurais mené une analyse univariée au travers d'un filtre différentiel pour retenir uniquement les variables qui sont le plus liées à ta variable réponse Behavior puis j'aurais appliqué ma régression logistique ou carrément une analyse discriminante de Fisher en fonction des hypothèses de normalité.

Ensuite, concernant la régression logistique polytomique (alors là il s'agit de souvenir mais ce week-end je peux prendre le temps de relire mes vieilles notes) elle fonctionne comme une régression logistique standard mais avec un système de pivot, cad, on fixe une classe de référence et on étudie chaque probabilité de classer dans une autre classe par rapport à elle.
Mais cette étude se fait une classe par une classe et non toutes les autres classes par rapport à la classe de référence. Il s'ensuit que la prédiction se fait en fonction de la probabilité la plus haute pour chacune d'elle. En fait de mémoire, car il y a 2 approches pour cette méthodo, vis à vis du package VGAM, la fonction va calculer les différentes probas en fonction des différentes classes par rapport au pivot et dés que ta proba dépasse 0.5 alors tu as la classe de prédiction, ex:

P[classe 2 | pivot= classe 1 / X] = 0.2
P[classe 3 | pivot= classe1 / X] = 0.6
P[classe 4 | pivot= classe 1 / X] = 0.7

Ici tu classeras ton observation en 3. (à confirmer mais il me semble que c'est un truc dans ce genre, le souci avec ce type de façon de faire c'est qu'on se base pas sur de la théorie, je parle bien de la règle de décision et non de l'élaboration des coefficients qui elle est fait sur des bases excellemment solides, mais sur du bidouillage).

Après, attention car on compte 2 types de régression logistique dans le cas de variable réponse à plus de 2 classes, celle qui suppose qu'il y a relation d'ordre et celle qui suppose que ce n'est pas le cas.

Bonjour et merci pour votre réponse!

Je vais essayer de suivre votre conseil. Cependant sauriez-vous me conseiller sur quelle package utiliser sous R?

J'avais bien compris les résultats de la régression, le soucis c'est que je ne sais pas quel comportement il prend comme référence et je ne sais pas comment le savoir avec ce package malheureusement.

Je vous siuhaites de très bonne fêtes!

Cordialement

Olivier Darbellay

Au vue de la sortie R que vous avez affiché je répondrais qu'il s'agit de la 4ème modalité de la variable réponse, surement celle la plus éloigné en suivant un ordre alphabétique.

Concernant le package, personnellement j'avais utilisé le package VGAM qui était assez simple d'utilisation et assez complet de souvenir.

C'est celui que j'ai utilisé, il y a la possibilité de faire une régression logistique simple avec ce package du coup je suppose ?

Faire une régression logistique binaire tu veux dire? J'avoue ne pas savoir... je pense que oui mais si ce n'est pas le cas tu peux utiliser la fonction glm qui est très complète (par contre prendre le temps de voir les fonctions annexes à utiliser sur un objet sortie de glm() pour pouvoir profiter réellement de tous les bien faits de cet outil) qui doit déjà être présente dans R.

Je ne vais pas pouvoir utiliser la régression binaire vu que ma variable a 4 états differents. Je vais utiliser la glm dans ce cas.
Milles merci et bonnnes fetes a vous!

Bonjour,

Deux choses: c'est très agréable de lire un premier message aussi bien présenté que le tien et ça donne envie de lire la suite.
Et pour le package se concentrer (de mémoire) sur le glm() qui est très complet et bien renseigné sur le net.

Bon courage pour la suite!

Kali

De souvenir glm() ne fait pas les régressions polytomiques mais uniquement les régressions logistiques binaires.

salut,

je confirme, la fonction glm ne propose pas de modèle multinomial via son argument family. Il ya beaucoup d'autres packages qui font ça comme par exemple le package nnet ou encore mlogit. Le package VGAM n'est pas tout à fait construit pour ça même si les fonctions à l'intérieur peuvent parfaitement convenir.

Par contre dans R la modalité de référence est par défaut toujours la première modalité dans l'ordre alphanumérique. Pour toi, le comportement de référence va donc être "Feeding".
Il est possible de modifier ça via les fonction contr.treatment() ou encore C() ou encore contrasts(). Pour cette dernière il faut définir une matrice de contrastes en mettant des 0 sur toute la ligne de la modalité de référence. Tu peux regarder la matrice des contrastes de ton facteur de comportement en faisant :

Code:

contrasts(Behavior)

nik

Merci pour les précisions, je pensais de mémoire qu'avec le package glm() c'était possible.

Kali

Merci à tous pour vos précisions, je vais continuer sur ma lancée alors, merci pour tout

lapive

Bonjour,

J'ai modifié mes données pour utiliser le package mlogit,

malheureusement, j'arrive bien a avoir mes données dans le format désiré

Code:

mlogit <- mlogit.data(Merge, choice = "Choice", shape = "long", alt.var = "Comp", drop.index = TRUE)

par contre quand je veux faire ma régression j'ai une erreure et je ne sais pas ce qu'elle veux dire, quelqu'un aurait une idée ?

Code:

> results <- mlogit(Choice ~ Temp + Valley + Age , mlogit)
Erreur dans solve.default(H, g[!fixed]) :
  routine Lapack dgesv : le système est exactement singulier : U[4,4] = 0
De plus : Il y a eu 12 avis (utilisez warnings() pour les visionner)

ça veut dire que la matrice du modèle n'est pas inversible car elle comporte au moins une colonne qui est une combinaison linéaire de 2 ou plusieurs autres colonnes.
Donc ici c'est le plan d'expérience qui est peut être mal conçu et qui est n'est pas équilibré.

nik

Ok,

J'ai fait un changement dans la formule (le changement suivant) :

Code:

results <- mlogit(Choice ~ 1 |Temp + Valley + Age + kid + Year + Month , mlogit)


et ca me donne un résultat :

Code:

> summary(results)

Call:
mlogit(formula = Choice ~ 1 | Temp + Valley + Age + kid + Year +
    Month, data = mlogit, method = "nr", print.level = 0)

Frequencies of alternatives:
        F        M        R      RUN
0.5011492 0.1144903 0.3778441 0.0065165

nr method
9 iterations, 0h:1m:42s
g'(-H)^-1g = 4.09E-05
successive function values within tolerance limits

Coefficients :
                      Estimate  Std. Error  t-value  Pr(>|t|)   
M:(intercept)      2.7522e+02  7.5105e+00  36.6441 < 2.2e-16 ***
R:(intercept)      8.3219e+01  5.5711e+00  14.9375 < 2.2e-16 ***
RUN:(intercept)    -1.3006e+02  5.0138e+01  -2.5941 0.0094844 **
M:Temp              2.0356e-02  9.5130e-04  21.3980 < 2.2e-16 ***
R:Temp            -4.9562e-02  7.0708e-04 -70.0939 < 2.2e-16 ***
RUN:Temp          -3.8838e-02  4.4967e-03  -8.6369 < 2.2e-16 ***
M:ValleyTrupchun    1.9034e-03  1.6126e-02  0.1180 0.9060407   
R:ValleyTrupchun    3.6090e-02  1.0385e-02  3.4752 0.0005105 ***
RUN:ValleyTrupchun  1.2421e+00  5.8255e-02  21.3215 < 2.2e-16 ***
M:Age              2.3546e-02  2.6771e-03  8.7955 < 2.2e-16 ***
R:Age              2.7134e-02  1.6988e-03  15.9721 < 2.2e-16 ***
RUN:Age            6.3920e-02  1.1322e-02  5.6455 1.647e-08 ***
M:kidY            -1.0488e-01  2.0389e-02  -5.1441 2.689e-07 ***
R:kidY            -4.7724e-01  1.2727e-02 -37.4974 < 2.2e-16 ***
RUN:kidY            5.7597e-01  6.9494e-02  8.2880 2.220e-16 ***
M:Year            -1.3769e-01  3.7381e-03 -36.8348 < 2.2e-16 ***
R:Year            -4.1159e-02  2.7718e-03 -14.8491 < 2.2e-16 ***
RUN:Year            5.9138e-02  2.4940e-02  2.3712 0.0177313 * 
M:Month            -4.4177e-02  8.0378e-03  -5.4961 3.883e-08 ***
R:Month            5.0182e-02  5.3523e-03  9.3759 < 2.2e-16 ***
RUN:Month          9.1490e-01  3.6274e-02  25.2221 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -264700
McFadden R^2:  0.021594
Likelihood ratio test : chisq = 11684 (p.value = < 2.22e-16)


es-ce que cela veux dire quelque chose ? j'ai un peu de peine a interpréter les résultats.

Ce que j'ai fait comme changement dans le fichier de base, c'est que j'ai mis des "yes" et "no" à la place de Y et N dans la colonne "choix" pour qu'il reconnaisse et mette "false et true" à la place de Y et N. Je ne sais pas si vous voyez ce que je veux dire mais ca m'a l'air de marcher.

Ce que j'ai fait comme changement dans le fichier de base, c'est que j'ai mis des "yes" et "no" à la place de Y et N dans la colonne "choix" pour qu'il reconnaisse et mette "false et true" à la place de Y et N

ça ne change rien. Tu peux lui mettre ce que tu veux comme nom de catégorie, s'il y en a 2 alors R codera toujours de la même façon.


Ce que tu présentes, ce sont les tests des paramètres pour savoir s'ils sont significativement différents de 0. ça n'a donc que peut d'intérêt par rapport à la question de savoir si une variable explicative a un effet ou non.

Il faut regarder s'il existe une méthode anova pour mlogit et c'est ça qu'il faut utiliser en construisant des modèles emboités et en les testant les uns par rapports aux autres. Normalement, ça devrait réaliser un test du rapport des vraisemblances.
Ce que tu présentes, donne par contre la valeur des paramètres c'est à dire le paramètre pour chaque variable explicative en fonction de chaque catégorie à expliquée toujours par comparaison à la catégorie témoin (feeding).
Par exemple, pour M:age, lorsque tu augmente l'âge d'une unité, le log du rapport des chances (p/1-p) entre la catégorie M et la catégorie Feeding augmente de 0.023, soit, exp(0.023)=1.023, 2.3% d'augementation.

Nik

Encore une fois merci pour ta contribution!

Ce que je fais la c'est que j'extrais les valeur des AIC pour chacun de mes modèles et je vais enlever des paramètres petit à petit pour maximiser la valeur, dès que j'ai la plus petite valeur d'AIC ca sera le modèle final que je vais utiliser dans mes résultats.

J'espère que ca va aller comme ca.

Sinon que pense-tu de la librairie lme4 avec comme random effect les individus ? après je ferais la même procédure, c'est à dire comparer les AIC entre les différent modèle

EDIT : je viens d'essayer avec lme4 (formule suivante :

Code:

fm1<-lmer(Behavior~1 + (1|Individual) + kid + Temp + Month + Age + Year, data=Mergeok)


)

et ce que je viens de trouver c'est un aic plus petit pour ca :

Code:

fm5<-lmer(Behavior~1 + (1|Individual) + Valley + kid + Temp + Month + Year, data=Mergeok)

alors que pour mlogit j'ai un AIC plus petit lorsque tous les paramètres sont présents dans la formule, quel est celui sur lequel je devrais me fier ?

d'autant que pour lme4, je n'ai pas essayé d'enlever d'autres paramètres pour la formule avec le plus petit AIC mais je ne serais pas étonné d'avoir des AIC plus petit lorsque j'enlève encore d'autres paramètres.

En tous cas merci beaucoup pour ton aide

je ne suis pas sûr que ton modèle soit bon dans lmer. tu ne spécifie pas la famille de distribution donc je ne sais pas s'il prend en compte le fait que behaviour est multinomial. Fais bien attention car chaque librairie à ses propres codes et la notation ne veut pas toujours dire la même chose. En plus il est fort probable que mlogit estime par MV alors que lmer estime via REML donc les AIC ne sont pas comparables.

Ta procédure de sélection de modèle n'est pas bonne car tu oublies le fait que l'ordre dans lequel tu vas éliminer les variables compte aussi dans la valeur de l'AIC. voir ?stepAIC

Pour résumer, tu es dans des librairies très complexes alors avant de partir dans toutes les directions, tu ferais mieux de te poser et de bien réfléchir aux questions que tu te poses sur tes données et ensuite de choisir les analyses adaptées et ne plus en sortir.

Nik a écrit:
Ta procédure de sélection de modèle n'est pas bonne car tu oublies le fait que l'ordre dans lequel tu vas éliminer les variables compte aussi dans la valeur de l'AIC. voir ?stepAIC

Je vais regarder ceci. Cependant j'ai enlevé chaque paramètres à la main.
Je vais aussi rester sur la librairie mlogit du coup. J'aurais bien voulu attendre un moment mais malheureusement le temps m'est compté et je n'ai pas eu d'aide de mon groupe concernant les statistiques, j'espère que ce que j'ai fait suffira à satisfaire mes superviseurs.

Je vous remercie pour votre aide précieuse en tous cas!

Je vais aussi rester sur la librairie mlogit du coup

c'est plus sage en effet.
utilise stepAIC, qui te fournir le "meilleur" modèle selon les données. Ne rentre pas n'importe quoi comme variable non plus sous peine de faire un modèle qui ne correspond pas à grand chose même s'il est le meilleur selon l'AIC


Nik

J'ai essayé la formule stepwise mais ca me sort une erreure a chaque fois :

Code:

> m1.step<-stepAIC(m1)
Start:  AIC=529131
Choice ~ 1 + (Temp + Valley + Age + Individual + kid + Year +
    Month)

Erreur dans solve.default(H, g[!fixed]) :
  routine Lapack dgesv : le système est exactement singulier : U[4,4] = 0


Donc je vais rester sur ce que j'ai fait à la main (j'ai enlevé 1 seul paramètre à la fois pour essayer toutes les combinaisons possibles)

[quote="Nik"]

Par exemple, pour M:age, lorsque tu augmente l'âge d'une unité, le log du rapport des chances (p/1-p) entre la catégorie M et la catégorie Feeding augmente de 0.023, soit, exp(0.023)=1.023, 2.3% d'augementation.

Nik

Salut, j'ai une dernière question

Je ne comprend pas comment tu es arrivé au 2.3% à partir de l'exp(0.023).
Faut-il simplement multiplier par 100 (pour avoir des %) ce que je trouve dans la colonne estimate ?

Ensuite, cela veux dire que si l'age augmente d'une unité j'ai 2.3% de chance en plus d'atteindre la catégorie M en partant de la catégorie Feeding ?

voila je veux juste m'assurer d'avoir bien compris pour avoir quelque chose qui tienne la route dans mon écrit.

Merci encore!

Lapive

exp(0.023) = 1.023, en % ça donne 102,3%. donc 2,3% d'augmentation par rapport à la proba de base. le%, c'est juste une façon de présenter les choses, il n'y a aucun calcul particulier à faire.

Ensuite, cela veux dire que si l'age augmente d'une unité j'ai 2.3% de chance en plus d'atteindre la catégorie M en partant de la catégorie Feeding ?

oui

Génial, merci pour tout, vous m'avez grandement facilité mon travail!

Lapive.

Edit :

J'ai encore une question, désolé, si la probabilité est négative (disons que j'ai -3.6% de passer d'un état a un autre lorsqu'un facteur augmente d'une unité) cela veux dire qu'il faut diminuer le facteur d'une unité pour avoir plus de chance de passer de l'état initial au nouvel état (Feeding -> moving par exemple)?
Et j'ai deux paramètres à deux états (Valley et kid) qu'il m'a affiché en tant que "ValleyTrupchun" (une des deux vallées) et "KidY" (l'autre cas c'est N pour no), étant donné qu'ils ne peuvent pas augmenter d'une unité, comment cela s'analyse ? c'est simplement le passage d'une vallée à l'autre ou la différence entre 1 enfant et pas d'enfant ?

Désolé de spam ce post mais je ne sais pas si vous aviez vous mon edit, malheureusement j'avais encore une question:

lapive a écrit:
J'ai encore une question, désolé, si la probabilité est négative (disons que j'ai -3.6% de passer d'un état a un autre lorsqu'un facteur augmente d'une unité) cela veux dire qu'il faut diminuer le facteur d'une unité pour avoir plus de chance de passer de l'état initial au nouvel état (Feeding -> moving par exemple)?
Et j'ai deux paramètres à deux états (Valley et kid) qu'il m'a affiché en tant que "ValleyTrupchun" (une des deux vallées) et "KidY" (l'autre cas c'est N pour no), étant donné qu'ils ne peuvent pas augmenter d'une unité, comment cela s'analyse ? c'est simplement le passage d'une vallée à l'autre ou la différence entre 1 enfant et pas d'enfant ?

Merci encore

Lapive