Przewidywanie uporządkowanego loginu w R.

Próbuję wykonać uporządkowaną regresję logit. Korzystam z takiego modelu (tylko głupi, mały model szacujący liczbę firm na rynku na podstawie miar dochodów i populacji). Moje pytanie dotyczy prognoz.

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

Kiedy uruchamiam przewidywanie (którego próbuję użyć, aby uzyskać przewidywaną wartość y), dane wyjściowe wynoszą 0, 3 lub 27, co w żaden sposób nie odzwierciedla prognozy opartej na moich manualnych prognozach ze współczynnika oszacowania i przechwyty. Czy ktoś wie, jak uzyskać „dokładne” prognozy dla mojego zamówionego modelu logit?

EDYTOWAĆ

Aby wyjaśnić moje obawy, moje dane odpowiedzi zawierają obserwacje na wszystkich poziomach

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

gdzie, jak moja przewidywana zmienna, wydaje się grupować

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

r econometrics logit ordered-logit

— prototoast
źródło

To jest dość niejasne. W jaki sposób wartości zwracane przez predictfunkcję różnią się od wartości wygenerowanych ręcznie? Jaka jest struktura twojej zmiennej zależnej? Podaj powtarzalny przykład.

— Sven Hohenstein

Myślę, że chciałbyś zobaczyć to-stats.stackexchange.com/questions/18119/…

— Blain Waan,

Nie do końca podążam za twoją sytuacją. Mówisz, że używasz modelu regresji porządkowej, ale także, jak najlepiej rozumiem, że zmienną odpowiedzi jest liczba firm na rynku. To jest liczba , to porządek, ale OLR nie jest właściwym sposobem na modelowanie tego; chcesz użyć jakiegoś wariantu regresji Poissona.

— gung - Przywróć Monikę

@gung Tak, rozumiem punkt o liczeniu vs. porządkowym. W tej chwili próbuję powielić papierowe pomysły.repec.org/ a/ ucp/jpolec/v99y1991i5p977-1009.html i używają regresji porządkowej. Oszacowałem również liczbę modeli, ale to nie pomaga mi w tym konkretnym zadaniu. Nie, nie chodzi o to, że po prostu chcę, aby R to zrobił, staram się zrozumieć, gdzie zachowanie odbiega od moich oczekiwań (ponieważ podejrzewam, że błąd jest z mojej strony, a nie R).

— prototoast

Czy zweryfikowałeś polr()względem innych funkcji? Można spróbować lrm()z pakietu rms: lrmFit <- lrm(y ~ pop0 + inc0); predict(lrmFit, type="fitted.ind"). Inną opcją jest vglm()z pakietu VGAM: vglmFit <- vglm(y ~ pop0 + inc0, family=propodds); predict(vglmFit, type="response"). Oba zwracają macierz przewidywanych prawdopodobieństw kategorii. Zobacz moją odpowiedź, aby uzyskać stamtąd przewidywane kategorie.

— caracal

Aby ręcznie zweryfikować przewidywania wynikające z użycia polr()z pakietu MASS, załóżmy sytuację z kategorycznie zależną zmienną z uporządkowanymi kategoriami oraz predyktorami . zakłada model proporcjonalnych szans $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})$

Aby zapoznać się z możliwymi wyborami zaimplementowanymi w innych funkcjach, zobacz tę odpowiedź . Funkcja logistyczna jest odwrotnością funkcji logit, więc przewidywane prawdopodobieństwa wynoszą $\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})}}$

Przewidywane prawdopodobieństwa kategorii to . Oto powtarzalny przykład w R z dwoma predyktorami . W przypadku porządkowej zmiennej przecięłem symulowaną zmienną ciągłą na 4 kategorie. $\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

Teraz dopasuj model proporcjonalnego prawdopodobieństwa za pomocą polr()i uzyskaj macierz przewidywanych prawdopodobieństw kategorii za pomocą 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

Aby ręcznie zweryfikować te wyniki, musimy wyodrębnić oszacowania parametrów, z nich obliczyć przewidywane logity, z tych logów obliczyć przewidywane prawdopodobieństwa , a następnie powiązać przewidywane prawdopodobieństwa kategorii z macierzą . $\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)

Porównaj z wynikiem z polr().

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

W przypadku przewidywanych kategorii predict(polr(), type="class")wystarczy wybrać - dla każdej obserwacji - kategorię o najwyższym prawdopodobieństwie.

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

Porównaj z wynikiem z 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
źródło