vcov.polr给出了负对角线元素(包括mwe)

Question

我对有序logits不太了解，所以也许我在做一些愚蠢的事情，但下面的内容似乎很简单，我不知道为什么vcov.polr会给我负对角线。我试着遵循vcov.polr中的代码，但是迷路了，所以我不知道出了什么问题。谢谢你的帮忙!

temp_fc <-
structure(list(the_index = c(4, 3, 3, 4, 4, 1, 3, 2, 3, 3, 3,
4, 4, 4, 1, 4, 4, 1, 4, 3, 4, 3, 3, 1, 4, 2, 4, 4, 4, 4, 3, 4,
1, 4, 4, 3, 4, 4, 2, 4, 1, 4, 4, 2, 3, 1, 4, 2, 1, 1, 4, 1, 4,
4, 4, 4, 1, 4, 3, 1, 4, 3, 4, 1, 4, 2, 3, 1, 4, 4, 1, 3, 2, 3,
2, 3, 3, 4, 1, 1, 4, 4, 4, 4, 1, 4, 4, 2, 4, 4, 4, 1, 1, 4, 3,
2, 4, 4, 1, 1), age = c(39L, 40L, 23L, 33L, 24L, 22L, 15L, 43L,
59L, 58L, 20L, 39L, 33L, 21L, 34L, 44L, 27L, 53L, 60L, 34L, 17L,
20L, 19L, 21L, 57L, 56L, 36L, 30L, 16L, 17L, 16L, 71L, 49L, 42L,
72L, 36L, 18L, 19L, 62L, 20L, 55L, 79L, 69L, 46L, 28L, 50L, 40L,
21L, 57L, 31L, 38L, 60L, 56L, 59L, 21L, 26L, 53L, 49L, 51L, 59L,
32L, 56L, 25L, 46L, 42L, 70L, 65L, 50L, 42L, 53L, 29L, 63L, 34L,
29L, 65L, 31L, 70L, 52L, 76L, 37L, 24L, 28L, 45L, 47L, 48L, 55L,
65L, 28L, 28L, 31L, 24L, 41L, 43L, 34L, 69L, 65L, 48L, 24L, 32L,
18L)), .Names = c("the_index", "age"), row.names = c("1006",
"1040", "1085", "1115", "1130", "1144", "1162", "1164", "1176",
"1178", "1181", "1192", "1223", "1258", "1259", "1307", "1311",
"1330", "1338", "1368", "1401", "1411", "1443", "1444", "1515",
"1538", "1600", "1628", "1632", "1683", "1692", "1701", "1710",
"1784", "1790", "1793", "1799", "1873", "1897", "1904", "1908",
"1982", "1988", "2015", "2029", "2052", "2065", "2077", "2088",
"2095", "2163", "10345", "10355", "10372", "10388", "10397",
"10414", "10415", "10417", "10421", "10428", "10430", "10456",
"10480", "10490", "10492", "10509", "10587", "10600", "10607",
"10609", "10617", "10625", "10626", "10630", "10640", "10674",
"10684", "10686", "10687", "10700", "10703", "10713", "10730",
"10740", "10747", "10750", "10762", "10792", "10794", "10803",
"10820", "10824", "10830", "10833", "10835", "10857", "10888",
"10917", "10933"), class = "data.frame")

library(MASS)
model_temp <-factor(the_index) ~ age + I(age^2)
the_OL <- polr(model_temp,temp_fc,method="logistic", Hess=TRUE)
diag(vcov(the_OL)) #this gives negative values
#so of course the following will have an error:
summary(the_OL)

Answer 1

我认为这是由于您在模型公式中使用了非正交多项式；这似乎导致了计算问题。使用你的代码，我得到了：

> the_OL2 <- polr(model_temp,temp_fc,method="logistic", Hess=TRUE)
> AIC(the_OL2)
[1] 254.6934

如果我通过poly()使用正交多项式重新拟合，我似乎得到了同样的拟合，但没有负方差的数值问题：

> the_OL <- polr(factor(the_index) ~ poly(age, 2), temp_fc, method="logistic", 
+  Hess=TRUE)
> AIC(the_OL)
[1] 254.6934
> diag(vcov(the_OL))
poly(age, 2)1 poly(age, 2)2           1|2           2|3 
   3.38608675    3.65981581    0.05958267    0.04704916 
          3|4 
   0.04103139

报告的两个模型的残差是相同的，两个模型的拟合值几乎相同(接近实际情况)。