Vorhersage des geordneten Logits in R.

Ich versuche eine geordnete Logit-Regression durchzuführen. Ich führe das Modell so aus (nur ein dummes kleines Modell, das die Anzahl der Unternehmen auf einem Markt anhand von Einkommens- und Bevölkerungsmaßen schätzt). Meine Frage betrifft Vorhersagen.

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

Wenn ich "Vorhersagen" ausführe (mit denen ich versuche, das vorhergesagte y zu erhalten), sind die Ausgaben entweder 0, 3 oder 27, was in keiner Weise die Vorhersage widerspiegelt, die auf meinen manuellen Vorhersagen aus dem Koeffizienten basieren sollte Schätzungen und Abschnitte. Weiß jemand, wie man "genaue" Vorhersagen für mein bestelltes Logit-Modell erhält?

BEARBEITEN

Um meine Bedenken zu verdeutlichen, enthalten meine Antwortdaten Beobachtungen auf allen Ebenen

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

wo wie meine Vorhersagevariable sich zu bündeln scheint

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

r econometrics logit ordered-logit Prototoast
quelle

Das ist ziemlich vage. Wie unterscheiden sich die von der predictFunktion zurückgegebenen Werte von den manuell generierten? Wie ist Ihre abhängige Variable aufgebaut? Bitte geben Sie ein reproduzierbares Beispiel an.

Sven Hohenstein

Ich denke, Sie möchten this-stats.stackexchange.com/questions/18119/…

Blain Waan

Ich verfolge deine Situation nicht ganz. Sie sagen, dass Sie ein ordinales Regressionsmodell verwenden, aber Sie sagen auch, wie ich am besten verstehe, dass Ihre Antwortvariable die Anzahl der Unternehmen in einem Markt ist. Das ist eine Zählung , es ist eine Ordnungszahl, aber OLR ist nicht der richtige Weg, dies zu modellieren. Sie möchten eine Variante der Poisson-Regression verwenden.

Gung - Reinstate Monica

@gung Ja, ich verstehe den Punkt über Zählung gegen Ordnungszahl. Im Moment versuche ich, das Papier ideas.repec.org/a/ucp/jpolec/v99y1991i5p977-1009.html zu replizieren, und sie verwenden eine ordinale Regression. Ich habe auch Zählmodelle geschätzt, aber das hilft mir bei dieser speziellen Aufgabe nicht. Nein, es ist nicht so, dass ich nur möchte, dass R dies tut. Ich versuche zu verstehen, wo das Verhalten von meinen Erwartungen abweicht (weil ich vermute, dass der Fehler meinerseits ist, nicht R).

Prototoast

Haben Sie sich polr()gegen andere Funktionen verifiziert ? Sie könnten lrm()aus Paket versuchen rms: lrmFit <- lrm(y ~ pop0 + inc0); predict(lrmFit, type="fitted.ind"). Eine weitere Option ist vglm()aus dem Paket VGAM: vglmFit <- vglm(y ~ pop0 + inc0, family=propodds); predict(vglmFit, type="response"). Beide geben die Matrix der vorhergesagten Kategoriewahrscheinlichkeiten zurück. Sehen Sie meine Antwort , um die vorhergesagten Kategorien von dort zu erhalten.

Caracal

Antworten:

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})$

$\hat{p}(Y \leqslant g)$

\hat{p} (Y. ⩽ 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})}}$

$\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

Passen Sie nun das Proportional-Odds-Modell mit an polr()und erhalten Sie die Matrix der vorhergesagten Kategoriewahrscheinlichkeiten mit 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

$\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)

Vergleiche mit dem Ergebnis von polr().

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

Wählen Sie für die vorhergesagten Kategorien predict(polr(), type="class")einfach - für jede Beobachtung - die Kategorie mit der höchsten Wahrscheinlichkeit aus.

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

Vergleichen Sie mit dem Ergebnis von polr().

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

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

Karakal
quelle