Comment comprendre la sortie de la fonction polr de R (régression logistique ordonnée)?

Je suis nouveau dans R, j'ai ordonné une régression logistique et polr.

La section "Exemples" au bas de la page d'aide de polr (qui adapte un modèle de régression logistique ou probit à une réponse factorielle ordonnée) montre

options(contrasts = c("contr.treatment", "contr.poly"))
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
pr <- profile(house.plr)
plot(pr)
pairs(pr)

Quelles informations prcontiennent-elles? La page d'aide sur le profil est générique et ne donne aucune indication pour polr.
Que plot(pr)montre- t-on ? Je vois six graphiques. Chacun a un axe X qui est numérique, bien que l'étiquette soit une variable indicatrice (ressemble à une variable d'entrée qui est un indicateur d'une valeur ordinale). Ensuite, l'axe Y est "tau", ce qui est complètement inexpliqué.
Que pairs(pr)montre- t-on ? Cela ressemble à un tracé pour chaque paire de variables d'entrée, mais là encore je ne vois aucune explication des axes X ou Y.
Comment comprendre si le modèle correspondait bien? summary(house.plr)montre la déviance résiduelle 3479.149 et AIC (Akaike Information Criterion?) de 3495.149. Est-ce bon? Dans le cas où celles-ci ne sont utiles que comme mesures relatives (c'est-à-dire pour les comparer à un autre modèle), qu'est-ce qu'une bonne mesure absolue? La déviance résiduelle est-elle approximativement distribuée en khi carré? Peut-on utiliser "% correctement prévu" sur les données d'origine ou une validation croisée? Quelle est la façon la plus simple de procéder?
Comment appliquer et interpréter anovasur ce modèle? Les documents disent "Il existe des méthodes pour les fonctions d'ajustement de modèle standard, y compris prédire, résumé, vcov, anova." Cependant, l'exécution anova(house.plr)entraîneanova is not implemented for a single "polr" object
Comment interpréter les valeurs de t pour chaque coefficient? Contrairement à certains ajustements de modèle, il n'y a pas de valeurs P ici.

Je me rends compte que c'est beaucoup de questions, mais il est logique pour moi de poser en un seul paquet ("comment puis-je utiliser cette chose?") Plutôt que 7 questions différentes. Toute information appréciée.

r logistic dfrankow
la source

@dfrankow Aide un peu grossière et certainement très partielle pour vos deux premières questions, mais methods("profile")vous donnera les méthodes (S3 dans ce cas) associées à un profileobjet R , alors vous verrez qu'il existe une méthode dédiée aux polrrésultats, que vous pouvez parcourir en ligne en tapant getAnywhere("profile.polr")à l'invite R.

chl

Merci! Le code source est bon. L'explication serait encore meilleure. :)

dfrankow

Quelqu'un m'a montré "Modern Applied Statistics with S" de Venables et Ripley. La section 7.3 contient «Un exemple de tableau de fréquences à quatre voies» qui couvre largement ce modèle de maison. Lecture ..

dfrankow

En fait, la section est "un modèle de cotes proportionnelles"

dfrankow

Réponses:

Je vous suggère de consulter des livres sur l'analyse des données catégorielles (cf. Analyse des données catégoriques d'Alan Agresti, 2002) pour une meilleure explication et compréhension de la régression logistique ordonnée . Toutes les questions que vous posez trouvent essentiellement leur réponse dans quelques chapitres de ces livres. Si vous n'êtes intéressé que par Rdes exemples connexes, Extending Linear Models in R par Julian Faraway (CRC Press, 2008) est une excellente référence.

Avant de répondre à vos questions, la régression logistique ordonnée est un cas de modèles logit multinomiaux dans lesquels les catégories sont ordonnées. Supposons que nous ayons commandé catégories et que , pour individu , en réponse ordinal , pour . Avec une réponse ordonnée, il est souvent plus facile de travailler avec les probabilités cumulées, $J$ $i$ $Y_i$ $p_{ij}=P(Yi=j)$ $j=1,..., J$ . Les probabilités cumulatives sont en augmentation et invariantes à la combinaison de catégories adjacentes. De plus, , nous n'avons donc besoin que desprobabilités dumodèle . $\gamma_{ij}=P(Y_i \le j)$ $\gamma_{iJ}=1$ $J–1$

Maintenant, nous voulons lier s aux covariables . Dans votre cas, a 3 niveaux ordonnés: , , . Il est plus logique de les traiter comme ordonnés plutôt que non ordonnés. Les variables restantes sont vos covariables. Le modèle spécifique que vous envisagez est le modèle des cotes proportionnelles et est mathématiquement équivalent à: $\gamma_{ij}$ $x$ Satlowmediumhigh

logit γ_{j} (x_{i}) = θ_{j} - β^{T} x_{i}, j = 1 \dots J - 1

$\mbox{logit } \gamma_j(x_i) = \theta_j - \beta^T x_i, j = 1 \ldots J-1$

where γ_{j} (x_{i}) = P (Y_{i} \leq j | x_{i})

$\mbox{where }\gamma_j(x_i)=P(Y_i \le j | x_i)$

On l'appelle ainsi parce que les cotes relatives pour comparant et sont: $Y \le j$ $x_1$ $x_2$

(\frac{γ_{j} (x_{1})}{1 - γ_{j} (x_{1})}) / (\frac{γ_{j} (x_{2})}{1 - γ_{j} (x_{2})}) = \exp (- β^{T} (x_{1} - x_{2}))

$\left(\frac {\gamma_j(x_1)}{1-\gamma_j(x_1)}\right) / \left(\frac {\gamma_j(x_2)}{1-\gamma_j(x_2)}\right)=\exp(-\beta^T (x_1-x_2))$

$j$

Maintenant, je vais répondre à quelques (1, 2, 4) questions.

Comment comprendre si le modèle correspondait bien? résumé (house.plr) montre la déviance résiduelle 3479.149 et AIC (Akaike Information Criterion?) de 3495.149. Est-ce bon? Dans le cas où celles-ci ne sont utiles que comme mesures relatives (c'est-à-dire pour les comparer à un autre modèle), qu'est-ce qu'une bonne mesure absolue? La déviance résiduelle est-elle approximativement distribuée en khi carré? Peut-on utiliser "% correctement prévu" sur les données d'origine ou une validation croisée? Quelle est la façon la plus simple de procéder?

Un modèle ajusté par polrest spécial glm, donc toutes les hypothèses qui valent pour une glmtenue traditionnelle ici. Si vous prenez bien soin des paramètres, vous pouvez déterminer la distribution. Plus précisément, pour tester si le modèle est bon ou non, vous pouvez faire un test d'adéquation , qui teste le null suivant (notez que c'est subtil, vous voulez surtout rejeter le null, mais ici vous ne voulez pas le rejeter pour obtenir un bon ajustement):

H_{o} : current model is good enough

$H_o: \mbox{ current model is good enough }$

Vous utiliseriez le test du chi carré pour cela. La valeur de p est obtenue comme:

1-pchisq(deviance(house.plr),df.residual(house.plr))

La plupart du temps, vous espérez obtenir une valeur de p supérieure à 0,05 afin de ne pas rejeter la valeur nulle pour conclure que le modèle est bon (l'exactitude philosophique est ignorée ici).

L'AIC doit être élevé pour un bon ajustement en même temps que vous ne voulez pas avoir un grand nombre de paramètres. stepAICest un bon moyen de vérifier cela.

Oui, vous pouvez certainement utiliser la validation croisée pour voir si les prédictions sont valables. Voir predictfonction (option:) type = "probs"dans ?polr. Il vous suffit de prendre soin des covariables.

Quelles informations contient pr? La page d'aide sur le profil est générique et ne donne aucune indication pour polr

Comme indiqué par @chl et autres, prcontient toutes les informations nécessaires pour obtenir des CI et d'autres informations liées à la probabilité du polr fit. Tous les glms sont ajustés en utilisant la méthode d'estimation des moindres carrés pondérée de manière itérative pour la vraisemblance logarithmique. Dans cette optimisation, vous obtenez beaucoup d'informations (veuillez consulter les références) qui seront nécessaires pour calculer la matrice de covariance de la variance, l'IC, la valeur t, etc. Elle comprend tout cela.

Comment interpréter les valeurs de t pour chaque coefficient? Contrairement à certains modèles> ajustements, il n'y a pas de valeurs P ici.

Contrairement au modèle linéaire normal (spécial glm), les autres glms n'ont pas la bonne distribution t pour les coefficients de régression. Par conséquent, tout ce que vous pouvez obtenir est les estimations des paramètres et leur matrice de covariance de variance asymptotique en utilisant la théorie du maximum de vraisemblance. Par conséquent:

Variance (\hat{β}) = (X^{T} W X)^{- 1} \hat{ϕ}

$\text{Variance}(\hat \beta) = (X^T W X)^{-1}\hat \phi$

L'estimation divisée par son erreur standard est ce que BDR et WV appellent la valeur t (je suppose la MASSconvention ici). Elle équivaut à la valeur t d'une régression linéaire normale mais ne suit pas une distribution t. En utilisant CLT, il est distribué normalement asymptotiquement. Mais ils préfèrent ne pas utiliser cette valeur approximative (je suppose), donc pas de valeurs p. (J'espère que je ne me trompe pas, et si je le suis, j'espère que BDR n'est pas sur ce forum. J'espère en outre que quelqu'un me corrigera si je me trompe.)

suncoolsu
la source

J'ajouterai plus.

suncoolsu

Merci pour cela. Je l'ai lu plusieurs fois. De nombreuses questions demeurent. 1. Fonctionnellement dans R, comment tester l'hypothèse de cotes proportionnelles? 2. Êtes-vous sûr que le test du chi carré est correct? Sur cet exemple, il renvoie 0, ce qui signifie .. ajustement merdique? Mais certaines des valeurs t sont assez élevées (InflHigh 10.1, InflMedium 5.4, ContHigh 3.7). 3. Que montrent les parcelles ou les paires?

dfrankow

Merci pour votre réponse détaillée suncoolsu. Je suis dans une situation similaire et j'ai quelques questions. 1. J'obtiens également 0 pour chaque modèle en utilisant votre équation de test chi carré. 2. La page wikipedia sur AIC dit "le modèle préféré est celui avec la valeur AIC minimale", mais vous avez dit, "AIC devrait être élevé pour un bon ajustement." J'essaie de concilier ces comptes.

Sam Swift

@dfrankow et @Sam Swift. Je suis désolé, j'ai été un peu occupé à écrire des papiers. Ok - si vous obtenez une valeur de p = 0, cela signifie que le modèle n'est PAS un bon ajustement car le test de qualité de l'ajustement échoue. Concernant le problème de l'AIC, wikipedia et moi utilisons une convention différente pour cela. J'utilise celui qui est utilisé par BDR et WV. (cf. Extension des modèles linéaires en R, par le Dr Julian Faraway)

suncoolsu

Il y a quelques questions dédiées pour les valeurs 0/1 p et l'interprétation AIC que vous pourriez trouver utiles: stats.stackexchange.com/questions/15223/… stats.stackexchange.com/questions/81427/…

Scott

J'ai beaucoup apprécié la conversation ici, mais je pense que les réponses n'ont pas correctement adressé tous les (très bons) éléments à la question que vous avez posée. La deuxième moitié de l'exemple de page polrconcerne le profilage. Une bonne référence technique ici est Venerables et Ripley qui discutent du profilage et de ce qu'il fait. Il s'agit d'une technique critique lorsque vous sortez de la zone de confort de l'ajustement des modèles de famille exponentielle avec la plus grande probabilité (GLM réguliers).

$k-1$ $k$ lmernls,, polret glm.nb.

La page d'aide de ?profile.glmdevrait être utile, car les polrobjets sont essentiellement des GLM (plus les seuils catégoriels). Enfin, vous pouvez réellement atteindre le code source, s'il est utile, en utilisant getS3method('profile', 'polr'). J'utilise getS3methodbeaucoup cette fonction car, alors que R semble insister sur le fait que de nombreuses méthodes doivent être cachées, on peut en apprendre beaucoup sur la mise en œuvre et les méthodes en examinant le code.

• Quelles informations contient pr? La page d'aide sur le profil est générique et ne donne aucune indication pour polr.

prest un profile.polr, profileobjet (classe héritée profile). Il y a une entrée pour chaque covariable. Le profileur fait une boucle sur chaque covariable, recalcule l'ajustement optimal du modèle avec cette covariable fixée à un montant légèrement différent. Le résultat montre la valeur fixe de la covariable mesurée sous la forme d'une différence "z-score" à l'échelle de sa valeur estimée et les effets fixes résultants dans d'autres covariables. Par exemple, si vous regardez pr$InflMedium, vous remarquerez que, lorsque "z" vaut 0, les autres effets fixes sont les mêmes que ceux trouvés dans l'ajustement d'origine.

• Que montre l'intrigue (pr)? Je vois six graphiques. Chacun a un axe X qui est numérique, bien que l'étiquette soit une variable indicatrice (ressemble à une variable d'entrée qui est un indicateur d'une valeur ordinale). Ensuite, l'axe Y est "tau", ce qui est complètement inexpliqué.

Encore une fois, ?plot.profiledonne la description. Le graphique montre approximativement comment les coefficients de régression covarient. tau est la différence mise à l'échelle, le score z avant, donc sa valeur 0 donne les coefficients d'ajustement optimaux, représentés avec une coche. Vous ne diriez pas que cet ajustement est si bien comporté, mais ces "lignes" sont en fait des splines. Si la probabilité se comportait de manière très irrégulière à l'ajustement optimal, vous observeriez un comportement étrange et imprévisible dans le graphique. Cela vous incomberait d'estimer la sortie à l'aide d'une estimation d'erreur plus robuste (bootstrap / jackknife), de calculer des IC à l'aide de method='profile', de recoder des variables ou d'effectuer d'autres diagnostics.

• Que montrent les paires (pr)? Cela ressemble à un tracé pour chaque paire de variables d'entrée, mais encore une fois, je ne vois aucune explication des axes X ou Y.

Le fichier d'aide dit: "La méthode des paires montre, pour chaque paire de paramètres x et y, deux courbes se coupant à l'estimation de vraisemblance maximale, qui donnent les loci des points auxquels les tangentes aux contours du profil bivarié deviennent verticales et horizontale, respectivement. Dans le cas d'une probabilité de profil normal exactement bivariée, ces deux courbes seraient des lignes droites donnant les moyennes conditionnelles de y | x et x | y, et les contours seraient exactement elliptiques. " Fondamentalement, ils vous aident à nouveau à visualiser les ellipses de confiance. Les axes non orthogonaux indiquent des mesures hautement covariables, telles que InfMedium et InfHigh qui sont intuitivement très liées. Encore une fois, des probabilités irrégulières conduiraient à des images assez déconcertantes ici.

• Comment peut-on comprendre si le modèle correspondait bien? résumé (house.plr) montre la déviance résiduelle 3479.149 et AIC (Akaike Information Criterion?) de 3495.149. Est-ce bon? Dans le cas où celles-ci ne sont utiles que comme mesures relatives (c'est-à-dire pour les comparer à un autre modèle), qu'est-ce qu'une bonne mesure absolue? La déviance résiduelle est-elle approximativement distribuée en khi carré? Peut-on utiliser "% correctement prévu" sur les données d'origine ou une validation croisée? Quelle est la façon la plus simple de procéder?

Une hypothèse qu'il est bon d'évaluer est l'hypothèse de cotes proportionnelles. Cela se reflète quelque peu dans le test global (qui évalue polr par rapport à un modèle log-linéaire saturé). Une limitation ici est qu'avec des données volumineuses, les tests globaux échouent toujours. Par conséquent, l'utilisation de graphiques et d'inspections d'estimations (bêtas) et de précision (ES) pour le modèle log-linéaire et l'ajustement polr est une bonne idée. S'ils sont massivement en désaccord, quelque chose ne va peut-être pas.

Avec des résultats ordonnés, il est difficile de définir un pourcentage d'accord. Comment choisirez-vous un classificateur basé sur le modèle, et si vous le faites, comment discuterez-vous des performances médiocres d'un classificateur médiocre. modeest un mauvais choix. Si j'ai 10 logits de catégorie et que ma prédiction est toujours mais une seule catégorie, ce n'est peut-être pas une mauvaise chose. De plus, mon modèle peut correctement prédire une probabilité de 40% de réponse 0, mais également une probabilité de 20% de 8, 9, 10. Donc, si j'observe 9, est-ce bon ou mauvais? Si vous devez mesurer l'accord, utilisez un kappa pondéré, ou même MSE. Le modèle log-linéaire produira toujours le meilleur accord. Ce n'est pas ce que fait le POLR.

• Comment appliquer et interpréter anova sur ce modèle? Les documents disent "Il existe des méthodes pour les fonctions d'ajustement de modèle standard, y compris prédire, résumé, vcov, anova." Cependant, l'exécution d'anova (house.plr) entraîne que anova n'est pas implémenté pour un seul objet "polr"

Vous pouvez tester des modèles imbriqués avec waldtestet lrtestdans le lmtestpackage dans R. C'est équivalent à ANOVA. L'interprétation est exactement la même que pour les GLM.

• Comment interprète-t-on les valeurs de t pour chaque coefficient? Contrairement à certains ajustements de modèle, il n'y a pas de valeurs P ici.

Encore une fois, contrairement aux modèles linéaires, le modèle POLR est capable d'avoir des problèmes avec une probabilité irrégulière, donc l'inférence basée sur la Hesse peut être très instable. C'est analogue à l'ajustement de modèles mixtes, voir par exemple le fichier d'aide sur confint.merModpour le paquet lme4. Ici, les évaluations faites avec le profilage montrent que la covariance se comporte bien. Les programmeurs l'auraient fait par défaut, sauf que le profilage peut être très intensif sur le plan des calculs, et donc ils le laissent entre vos mains. Si vous devez voir l'inférence basée sur Wald, utilisez-la coeftest(house.plr)dans le lrtestpackage.

AdamO
la source

Pour «tester» (c.-à-d. Évaluer) l'hypothèse de cotes proportionnelles dans R, vous pouvez utiliser residuals.lrm () dans le package de conception de Frank Harrell Jr. Si vous tapez? Residuals.lrm, il y a un exemple rapide à reproduire de la façon dont Frank Harrell recommande d'évaluer l'hypothèse de cotes proportionnelles (c'est-à-dire visuellement, plutôt que par un test à bouton-poussoir). La conception évalue les régressions logistiques ordonnées à l'aide de lrm (), que vous pouvez remplacer par polr () de MASS.

Pour un exemple plus formel de la façon de tester visuellement l'hypothèse de cotes proportionnelles dans R, voir: Document: Modèles de régression à réponse ordinale en écologie Auteur (s): Antoine Guisan et Frank E. Harrell Source: Journal of Vegetation Science, Vol. 11, n ° 5 (oct., 2000), p. 617-626

mBrewster
la source

J'apprécie sincèrement votre réponse. Cependant, le but de StackExchange est de fournir des réponses, pas des références. Les statisticiens semblent particulièrement enclins à ce problème de référence. Y a-t-il des détails que vous pouvez ajouter sur la façon d'utiliser residuals.lrm? Par exemple, un exemple de commande et un exemple d'interprétation du graphique pour l'exemple house.plr?

dfrankow

Mise à jour à partir du site Web de l'auteur: "Le package Design est désormais obsolète. Les utilisateurs R doivent utiliser le package rms à la place". Mark, votre réponse m'a été très utile.

Tal Galili