Прогнозирование заказанного логита в R

12

Я пытаюсь сделать упорядоченную регрессию логита. Я управляю моделью вот так (просто глупая маленькая модель, оценивающая количество фирм на рынке по показателям дохода и населения). Мой вопрос о прогнозах.

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

Когда я запускаю прогнозирование (которое я пытаюсь использовать для получения прогнозируемого значения y), выходные значения равны 0, 3 или 27, что никоим образом не отражает то, что должно казаться прогнозом, основанным на моих ручных прогнозах из коэффициента оценки и перехватывает. Кто-нибудь знает, как получить "точные" прогнозы для моей заказанной модели логита?

РЕДАКТИРОВАТЬ

Чтобы прояснить мою озабоченность, мои данные ответов содержат наблюдения на всех уровнях

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

где, как моя прогнозируемая переменная, кажется, накапливается

> head(table(pr_out))
pr_out
0     1   2   3   4   5 
117   0   0 114   0   0 
prototoast
источник
2
Это довольно расплывчато. Чем значения, возвращаемые predictфункцией, отличаются от значений , сгенерированных вами вручную? Какова структура вашей зависимой переменной? Пожалуйста, приведите воспроизводимый пример.
Свен Хоэнштейн
1
Я думаю , что вы хотели бы видеть this- stats.stackexchange.com/questions/18119/...
Блейн Waan
2
Я не совсем понимаю вашу ситуацию. Вы говорите, что используете модель порядковой регрессии, но вы также говорите, насколько я понимаю, что вашей переменной отклика является количество фирм на рынке. Это счет , это порядковый номер, но OLR не является правильным способом моделирования этого; Вы хотите использовать какой-либо вариант регрессии Пуассона.
gung - Восстановить Монику
2
@ Gung Да, я понимаю смысл счета против порядкового номера. В данный момент я пытаюсь воспроизвести статью ideas.repec.org/a/ucp/jpolec/v99y1991i5p977-1009.html, и они используют порядковый регресс. Я также оценил количество моделей, но это не помогает мне с этой конкретной задачей. Кроме того, нет, я не хочу, чтобы R просто делал это, я пытаюсь понять, где поведение отклоняется от моих ожиданий (потому что я подозреваю, что ошибка с моей стороны, а не с R).
протостаст
1
Вы проверяли polr()по другим функциям? Вы можете попробовать lrm()из пакета rms: lrmFit <- lrm(y ~ pop0 + inc0); predict(lrmFit, type="fitted.ind"). Другой вариант vglm()из пакета VGAM: vglmFit <- vglm(y ~ pop0 + inc0, family=propodds); predict(vglmFit, type="response"). Оба возвращают матрицу предсказанных вероятностей категории. Смотрите мой ответ, чтобы получить прогнозируемые категории оттуда.
Каракал

Ответы:

23

polr()MASSY1,,g,,kX1,,Xj,,Xppolr()

logit(p(Yg))=lnp(Yg)p(Y>g)=β0g(β1X1++βpXp)

Для возможных вариантов, реализованных в других функциях, см. Этот ответ . Логистическая функция является обратной по отношению к логит-функции, так что предсказанной вероятности являютсяp^(Yg)

p^(Yg)=eβ^0g(β^1X1++β^pXp)1+eβ^0g(β^1X1++β^pXp)

Предсказанные вероятности категории: . Вот воспроизводимый пример в R с двумя предикторами . Для порядковой переменной я разрезал моделируемую непрерывную переменную на 4 категории.Х1,Х2YP^(Y=g)=P^(Yg)P^(Yg1)X1,X2Y

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

Теперь подгоните модель пропорциональных коэффициентов, используя polr()и получите матрицу предсказанных вероятностей категории, используя 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

Чтобы вручную проверить эти результаты, нам нужно извлечь оценки параметров из этих вычислений предсказанных логитов, из этих логитов вычислить предсказанные вероятности , а затем связать предсказанные вероятности категории с матрицей ,p^(Yg)

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)

Сравните с результатом из polr().

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

Для прогнозируемых категорий predict(polr(), type="class")просто выберите - для каждого наблюдения - категорию с наибольшей вероятностью.

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

Сравните с результатом polr().

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

> all.equal(factor(categHat), facHat, check.attributes=FALSE)  # manual verification
[1] TRUE
каракал
источник