I'm trying to calculate the variance-covariance matrix of a polr model using the Hessian matrix output from the function. This is from the example in the polr help file.
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess=TRUE)
house.plr
vcov(house.plr)
Yields:
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
When I take the inverse of the Hessian, I get the same matrix except for the last row and column, which corresponds to the Medium|High intercept.
(solve(house.plr$Hessian))
Yields:
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
If I then try to derive the standard errors of the coefficients using the Hessian-derived variance-covariance matrix sqrt(diag(solve(house.plr$Hessian))) , I end up getting a different estimate of the standard error for the Medium|High intercept and I would get from summary(house.plr). Can someone explain what is happening here?
EDIT: aosmith's comment below provides the simple answer to this question which is that vcov.polr does not simply return the inverse Hessian, but also involves the estimated thresholds. My confusion, the explication of which may be of use to others, stemmed from the fact that using vcov on an object returned from the ordered logistic function svyolr in the survey package returns a matrix that does NOT account for the estimated thresholds. It appears as though vcov.svyolr simply returns the inverse of the Hessian, which is not sufficient to derive the variance-covariance matrix in ordered logistic models (someone please correct me if I'm wrong on this). Assuming that it did got me into trouble when I tried to use it in further analysis. (I ultimately got around this by stealing code from summary.svyolr, which calculates the correct vcov to determine standard errors.)
To show what I mean, I continue the above code by reformatting the dataset above so that it can be used by svyolr, and then attempt to extract the vcov matrix using '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
This returns a matrix that can be contrasted to the one returned by vcov(house.plr) above, note the last row and last column:
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
