Вычислить дисперсионно-ковариационную матрицу, используя гессиан в polr
Я пытаюсь вычислить дисперсионно-ковариационную матрицу модели polr, используя матрицу Гессиана, выведенную из функции. Это из примера в файле справки polr.
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess=TRUE)
house.plr
vcov(house.plr)
Урожайность:
InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace ContHigh Low|Medium Medium|High
InflMedium 0.0109522057 0.0058388698 -0.0001186328 0.0003571197 0.0001194372 0.0005232715 0.005775204 0.006214155
InflHigh 0.0058388698 0.0161686862 -0.0011185535 -0.0005789764 -0.0002820811 0.0014699005 0.005527505 0.006480153
TypeApartment -0.0001186328 -0.0011185535 0.0142177040 0.0096212170 0.0096814809 -0.0013601564 0.008675546 0.008249148
TypeAtrium 0.0003571197 -0.0005789764 0.0096212170 0.0240787651 0.0097933321 -0.0019949511 0.008588096 0.008323855
TypeTerrace 0.0001194372 -0.0002820811 0.0096814809 0.0097933321 0.0229480160 -0.0020860814 0.008860681 0.008043844
ContHigh 0.0005232715 0.0014699005 -0.0013601564 -0.0019949511 -0.0020860814 0.0091270890 0.004438238 0.004703586
Low|Medium 0.0057752039 0.0055275049 0.0086755464 0.0085880959 0.0088606811 0.0044382384 0.015586835 0.014367825
Medium|High 0.0062141551 0.0064801526 0.0082491485 0.0083238553 0.0080438438 0.0047035861 0.014367825 0.015743208
Когда я беру инверсию гессиана, я получаю ту же матрицу, за исключением последней строки и столбца, которая соответствует пересечению Medium|High.
(solve(house.plr$Hessian))
Урожайность:
InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace ContHigh Low|Medium Medium|High
InflMedium 0.0109522057 0.0058388698 -0.0001186328 0.0003571197 0.0001194372 0.0005232715 0.005775204 0.0003698475
InflHigh 0.0058388698 0.0161686862 -0.0011185535 -0.0005789764 -0.0002820811 0.0014699005 0.005527505 0.0008026733
TypeApartment -0.0001186328 -0.0011185535 0.0142177040 0.0096212170 0.0096814809 -0.0013601564 0.008675546 -0.0003592705
TypeAtrium 0.0003571197 -0.0005789764 0.0096212170 0.0240787651 0.0097933321 -0.0019949511 0.008588096 -0.0002226414
TypeTerrace 0.0001194372 -0.0002820811 0.0096814809 0.0097933321 0.0229480160 -0.0020860814 0.008860681 -0.0006882434
ContHigh 0.0005232715 0.0014699005 -0.0013601564 -0.0019949511 -0.0020860814 0.0091270890 0.004438238 0.0002235743
Low|Medium 0.0057752039 0.0055275049 0.0086755464 0.0085880959 0.0088606811 0.0044382384 0.015586835 -0.0010271027
Medium|High 0.0003698475 0.0008026733 -0.0003592705 -0.0002226414 -0.0006882434 0.0002235743 -0.001027103 0.0018418271
Если я тогда попытаюсь вывести стандартные ошибки коэффициентов, используя полученную по гессиану матрицу дисперсии-ковариации sqrt(diag(solve(house.plr$Hessian))) Я получаю другую оценку стандартной ошибки для перехвата Medium|High и получаю summary(house.plr), Может кто-нибудь объяснить, что здесь происходит?
РЕДАКТИРОВАТЬ: комментарий Aosmith ниже дает простой ответ на этот вопрос, который заключается в том, что vcov.polr не просто возвращает обратный гессиан, но и включает в себя расчетные пороги. Моя путаница, объяснение которой может быть полезным для других, проистекает из того факта, что использование vcov на объекте, возвращенном из упорядоченной логистической функции svyolr в опросном пакете возвращается матрица, которая НЕ учитывает предполагаемые пороговые значения. Кажется, будто vcov.svyolr просто возвращает инверсию гессиана, что недостаточно для получения матрицы дисперсии-ковариации в упорядоченных логистических моделях (кто-то, пожалуйста, исправьте меня, если я ошибаюсь в этом). Предполагая, что это доставило мне неприятности, когда я попытался использовать его в дальнейшем анализе. (В конечном итоге я обошел это путем кражи кода из summary.svyolr, который вычисляет правильный vcov для определения стандартных ошибок.)
Чтобы показать, что я имею в виду, я продолжаю приведенный выше код, переформатируя указанный выше набор данных, чтобы он мог использоваться svyolrи затем попытайтесь извлечь матрицу vcov, используя 'vcov':
#Formatting
housing.long<-expand.table(housing, freq="Freq")
housing.long$Sat<-ordered(housing.long$Sat, levels=c("Low","Medium","High"))
housing.long$Infl<-relevel(housing.long$Infl,ref="Low")
housing.long$Type<-relevel(housing.long$Type,ref="Tower")
housing.long$Cont<-relevel(housing.long$Cont,ref="Low")
#create svy objecct
house.plrsvy <- svydesign(ids=~1, data=housing.long)
#store model
house.plrsvy <-svyolr(Sat ~ Infl + Type + Cont, design=house.plrsvy)
vcov(house.plrsvy) [c(2,1,3,4,5,6,7,8),c(2,1,3,4,5,6,7,8)] #reording so it matches the matrix above
Это возвращает матрицу, которая может быть противопоставлена матрице, возвращенной vcov(house.plr) выше, обратите внимание на последний ряд и последний столбец:
InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace ContHigh Low|Medium Medium|High
InflMedium 0.0107386871 0.0056386668 -0.0002981046 -0.0001050978 -0.0002336774 0.0004818157 0.005346583 0.0003733560
InflHigh 0.0056386668 0.0165482681 -0.0010957497 -0.0003494497 -0.0002874101 0.0012244473 0.005220762 0.0008554517
TypeApartment -0.0002981046 -0.0010957497 0.0146171435 0.0099883268 0.0100760460 -0.0011474964 0.009094801 -0.0003656346
TypeAtrium -0.0001050978 -0.0003494497 0.0099883268 0.0234161270 0.0102936420 -0.0020927259 0.008830475 -0.0003345245
TypeTerrace -0.0002336774 -0.0002874101 0.0100760460 0.0102936420 0.0231356022 -0.0023225443 0.008963416 -0.0006818656
ContHigh 0.0004818157 0.0012244473 -0.0011474964 -0.0020927259 -0.0023225443 0.0092707384 0.004737418 0.0001668623
Low|Medium 0.0053465831 0.0052207623 0.0090948014 0.0088304752 0.0089634162 0.0047374180 0.015934585 -0.0010762218
Medium|High 0.0003733560 0.0008554517 -0.0003656346 -0.0003345245 -0.0006818656 0.0001668623 -0.001076222 0.0018436991