vcov.polr is giving negative diagonal elements (mwe included)

I don't know too much about ordered logits so perhaps I'm doing something silly, but the following seems straightforward and I don't know why vcov.polr gives me negative diagonals. I tried to follow the code in vcov.polr but got lost so I don't know what is wrong. Thanks for your help!

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,开发者_Go百科 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)

---

I think this is due to your use of non-orthogonal polynomials in your model formulae; this appears to be causing the computational issues. With your code I get:

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

If I refit using orthogonal polynomials via poly(), I seem to get the same fit but without the numerical issues of negative variances:

> 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

The residual deviances reported for the two models are the same and the fitted values from the two models are almost the same (near as matters).