Prédire le logit ordonné dans R

J'essaie de faire une régression logit ordonnée. J'exécute le modèle comme ça (juste un petit modèle stupide qui estime le nombre d'entreprises sur un marché à partir des mesures du revenu et de la population). Ma question concerne les prédictions.

nfirm.opr<-polr(y~pop0+inc0, Hess = TRUE)
pr_out<-predict(nfirm.opr)

Lorsque je lance Predict (que j'essaie d'utiliser pour obtenir le y prévu), les sorties sont soit 0, 3 ou 27, ce qui ne reflète en rien ce qui devrait sembler être la prédiction basée sur mes prédictions manuelles du coefficient estimations et interceptions. Est-ce que quelqu'un sait comment obtenir des prévisions "précises" pour mon modèle logit commandé?

ÉDITER

Pour clarifier ma préoccupation, mes données de réponse contiennent des observations à tous les niveaux

>head(table(y))
y
0  1  2  3  4  5 
29 21 19 27 15 16

où comme ma variable prédite semble se regrouper

> head(table(pr_out))
pr_out
0     1   2   3   4   5 
117   0   0 114   0   0

r econometrics logit ordered-logit prototoast
la source

C'est assez vague. En quoi les valeurs renvoyées par la predictfonction diffèrent-elles de celles que vous avez générées manuellement? Quelle est la structure de votre variable dépendante? Veuillez fournir un exemple reproductible.

Sven Hohenstein

Je pense que vous voudriez voir ça- stats.stackexchange.com/questions/18119/…

Blain Waan

Je ne suis pas tout à fait à votre place. Vous dites que vous utilisez un modèle de régression ordinale, mais vous dites aussi, d'après ce que je comprends, que votre variable de réponse est le nombre d'entreprises sur un marché. C'est un décompte , c'est ordinale, mais OLR n'est pas la bonne façon de modéliser cela; vous souhaitez utiliser une variante de la régression de Poisson.

gung - Rétablir Monica

@gung Oui, je comprends le point sur le nombre vs ordinal. Pour le moment, j'essaie de reproduire le papier ideas.repec.org/a/ucp/jpolec/v99y1991i5p977-1009.html et ils utilisent une régression ordinale. J'ai également estimé des modèles de comptage, mais cela ne m'aide pas dans cette tâche particulière. De plus, non, ce n'est pas que je veux juste que R fasse ça, j'essaie de comprendre où le comportement s'écarte de mes attentes (parce que je soupçonne que l'erreur est de ma part, pas R).

prototoast

Avez-vous vérifié par polr()rapport à d'autres fonctions? Vous pouvez essayer lrm()de l' emballage rms: lrmFit <- lrm(y ~ pop0 + inc0); predict(lrmFit, type="fitted.ind"). Une autre option est vglm()de l' emballage VGAM: vglmFit <- vglm(y ~ pop0 + inc0, family=propodds); predict(vglmFit, type="response"). Les deux renvoient la matrice des probabilités de catégorie prédites. Voir ma réponse pour obtenir les catégories prévues à partir de là.

caracal

Réponses:

polr()MASS $Y$ $1, \ldots, g, \ldots, k$ $X_{1}, \ldots, X_{j}, \ldots, X_{p}$ polr()

logit (p (Y ⩽ g)) = \ln \frac{p (Y ⩽ g)}{p (Y > g)} = β_{0_{g}} - (β_{1} X_{1} + \dots + β_{p} X_{p})

$\text{logit}(p(Y \leqslant g)) = \ln \frac{p(Y \leqslant g)}{p(Y > g)} = \beta_{0_g} - (\beta_{1} X_{1} + \dots + \beta_{p} X_{p})$

Pour les choix possibles implémentés dans d'autres fonctions, voir cette réponse . La fonction logistique est l'inverse de la fonction logit, donc les probabilités prédites sont $\hat{p}(Y \leqslant g)$

\hat{p} (Oui ⩽ g) = \frac{e^{{\hat{β}}_{0_{g}} - ({\hat{β}}_{1} X_{1} + \dots + {\hat{β}}_{p} X_{p})}}{1 + e^{{\hat{β}}_{0_{g}} - ({\hat{β}}_{1} X_{1} + \dots + {\hat{β}}_{p} X_{p})}}

$\hat{p}(Y \leqslant g) = \frac{e^{\hat{\beta}_{0_{g}} - (\hat{\beta}_{1} X_{1} + \dots + \hat{\beta}_{p} X_{p})}}{1 + e^{\hat{\beta}_{0_{g}} - (\hat{\beta}_{1} X_{1} + \dots + \hat{\beta}_{p} X_{p})}}$

Les probabilités de catégorie prédites sont . Voici un exemple reproductible en R avec deux prédicteurs . Pour une variable ordinale , j'ai découpé une variable continue simulée en 4 catégories. $\hat{P}(Y=g) = \hat{P}(Y \leq g) - \hat{P}(Y \leq g-1)$ $X_{1}, X_{2}$ $Y$

set.seed(1.234)
N     <- 100                                    # number of observations
X1    <- rnorm(N, 5, 7)                         # predictor 1
X2    <- rnorm(N, 0, 8)                         # predictor 2
Ycont <- 0.5*X1 - 0.3*X2 + 10 + rnorm(N, 0, 6)  # continuous dependent variable
Yord  <- cut(Ycont, breaks=quantile(Ycont), include.lowest=TRUE,
             labels=c("--", "-", "+", "++"), ordered=TRUE)    # ordered factor

Ajustez maintenant le modèle de cotes proportionnelles en utilisant polr()et obtenez la matrice des probabilités de catégorie prédites en utilisant predict(polr(), type="probs").

> library(MASS)                              # for polr()
> polrFit <- polr(Yord ~ X1 + X2)            # ordinal regression fit
> Phat    <- predict(polrFit, type="probs")  # predicted category probabilities
> head(Phat, n=3)
         --         -         +        ++
1 0.2088456 0.3134391 0.2976183 0.1800969
2 0.1967331 0.3068310 0.3050066 0.1914293
3 0.1938263 0.3051134 0.3067515 0.1943088

Pour vérifier manuellement ces résultats, nous devons extraire les estimations des paramètres, à partir de celles-ci calculer les logits prédits, à partir de ces logits calculer les probabilités prédites , puis lier les probabilités de catégorie prédites à une matrice . $\hat{p}(Y \leqslant g)$

ce <- polrFit$coefficients         # coefficients b1, b2
ic <- polrFit$zeta                 # intercepts b0.1, b0.2, b0.3
logit1 <- ic[1] - (ce[1]*X1 + ce[2]*X2)
logit2 <- ic[2] - (ce[1]*X1 + ce[2]*X2)
logit3 <- ic[3] - (ce[1]*X1 + ce[2]*X2)
pLeq1  <- 1 / (1 + exp(-logit1))   # p(Y <= 1)
pLeq2  <- 1 / (1 + exp(-logit2))   # p(Y <= 2)
pLeq3  <- 1 / (1 + exp(-logit3))   # p(Y <= 3)
pMat   <- cbind(p1=pLeq1, p2=pLeq2-pLeq1, p3=pLeq3-pLeq2, p4=1-pLeq3)  # matrix p(Y = g)

Comparez avec le résultat de polr().

> all.equal(pMat, Phat, check.attributes=FALSE)
[1] TRUE

Pour les catégories prédites, il predict(polr(), type="class")suffit de sélectionner - pour chaque observation - la catégorie ayant la probabilité la plus élevée.

> categHat <- levels(Yord)[max.col(Phat)]   # category with highest probability
> head(categHat)
[1] "-"  "-"  "+"  "++" "+"  "--"

Comparez pour résulter de polr().

> facHat <- predict(polrFit, type="class")  # predicted categories
> head(facHat)
[1] -  -  +  ++ +  --
Levels: -- - + ++

> all.equal(factor(categHat), facHat, check.attributes=FALSE)  # manual verification
[1] TRUE

caracal
la source