>  library(MASS)
>  data(housing)
>  names(housing)
[1] "Sat"  "Infl" "Type" "Cont" "Freq"
>  sapply(housing, levels)
$Sat
[1] "Low"    "Medium" "High"

$Infl
[1] "Low"    "Medium" "High"

$Type
[1] "Tower"     "Apartment" "Atrium"    "Terrace"

$Cont
[1] "Low"  "High"

$Freq
NULL

>
>  out0 <- glm(Freq ~ Infl * Type * Cont + Sat,
+      family = poisson, data = housing)
>
>  addterm(out0, ~ . + Sat : (Infl + Type + Cont), test = "Chisq")
Single term additions

Model:
Freq ~ Infl * Type * Cont + Sat
         Df Deviance    AIC    LRT   Pr(Chi)
<none>        217.46 610.43
Infl:Sat  4   111.08 512.05 106.37 < 2.2e-16 ***
Type:Sat  6   156.79 561.76  60.67 3.292e-11 ***
Cont:Sat  2   212.33 609.30   5.13   0.07708 .
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>  out1 <- update(out0, . ~ . + Sat : (Infl + Type + Cont))
>  summary(out1, cor = FALSE)

Call:
glm(formula = Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont +
    Type:Cont + Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont,
    family = poisson, data = housing)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.6022  -0.5282  -0.0641   0.5757   1.9322

Coefficients:
                                   Estimate Std. Error z value Pr(>|z|)
(Intercept)                        3.135074   0.120094  26.105  < 2e-16 ***
InflMedium                         0.248327   0.159949   1.553 0.120532
InflHigh                          -0.412644   0.184923  -2.231 0.025652 *
TypeApartment                      0.292524   0.157457   1.858 0.063197 .
TypeAtrium                        -0.792847   0.214397  -3.698 0.000217 ***
TypeTerrace                       -1.018074   0.221208  -4.602 4.18e-06 ***
ContHigh                          -0.001407   0.169668  -0.008 0.993384
Sat.L                             -0.098106   0.112561  -0.872 0.383438
Sat.Q                              0.285658   0.122272   2.336 0.019479 *
InflMedium:TypeApartment          -0.017882   0.210468  -0.085 0.932292
InflHigh:TypeApartment             0.386869   0.233275   1.658 0.097231 .
InflMedium:TypeAtrium             -0.360311   0.304954  -1.182 0.237395
InflHigh:TypeAtrium               -0.036788   0.334727  -0.110 0.912486
InflMedium:TypeTerrace             0.185154   0.288803   0.641 0.521451
InflHigh:TypeTerrace               0.310748   0.334759   0.928 0.353265
InflMedium:ContHigh               -0.200060   0.228704  -0.875 0.381706
InflHigh:ContHigh                 -0.725790   0.282316  -2.571 0.010145 *
TypeApartment:ContHigh             0.569691   0.212113   2.686 0.007236 **
TypeAtrium:ContHigh                0.702114   0.276025   2.544 0.010970 *
TypeTerrace:ContHigh               1.215929   0.269903   4.505 6.64e-06 ***
InflMedium:Sat.L                   0.519627   0.096813   5.367 7.99e-08 ***
InflHigh:Sat.L                     1.140302   0.118161   9.650  < 2e-16 ***
InflMedium:Sat.Q                  -0.064474   0.102660  -0.628 0.529983
InflHigh:Sat.Q                     0.115436   0.127791   0.903 0.366355
TypeApartment:Sat.L               -0.520170   0.109766  -4.739 2.15e-06 ***
TypeAtrium:Sat.L                  -0.288484   0.149523  -1.929 0.053686 .
TypeTerrace:Sat.L                 -0.998666   0.141467  -7.059 1.67e-12 ***
TypeApartment:Sat.Q                0.055418   0.118505   0.468 0.640043
TypeAtrium:Sat.Q                  -0.273820   0.149701  -1.829 0.067383 .
TypeTerrace:Sat.Q                 -0.032328   0.149230  -0.217 0.828496
ContHigh:Sat.L                     0.340703   0.087762   3.882 0.000104 ***
ContHigh:Sat.Q                    -0.097929   0.094062  -1.041 0.297822
InflMedium:TypeApartment:ContHigh  0.046901   0.286174   0.164 0.869819
InflHigh:TypeApartment:ContHigh    0.126229   0.338175   0.373 0.708951
InflMedium:TypeAtrium:ContHigh     0.157239   0.390686   0.402 0.687339
InflHigh:TypeAtrium:ContHigh       0.478611   0.444178   1.078 0.281248
InflMedium:TypeTerrace:ContHigh   -0.500162   0.367041  -1.363 0.172981
InflHigh:TypeTerrace:ContHigh     -0.463099   0.454652  -1.019 0.308403
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 833.657  on 71  degrees of freedom
Residual deviance:  38.662  on 34  degrees of freedom
AIC: 455.63

Number of Fisher Scoring iterations: 3

>  dropterm(out1, test = "Chisq")
Single term deletions

Model:
Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont +
    Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont
               Df Deviance    AIC    LRT   Pr(Chi)
<none>               38.66 455.63
Infl:Sat        4   147.78 556.75 109.12 < 2.2e-16 ***
Type:Sat        6   100.89 505.86  62.23 1.586e-11 ***
Cont:Sat        2    54.72 467.69  16.06 0.0003256 ***
Infl:Type:Cont  6    43.95 448.92   5.29 0.5072454
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> addterm(out1, ~ . + Sat : (Infl + Type + Cont)^2, test = "Chisq")
Single term additions

Model:
Freq ~ Infl + Type + Cont + Sat + Infl:Type + Infl:Cont + Type:Cont +
    Infl:Sat + Type:Sat + Cont:Sat + Infl:Type:Cont
              Df Deviance    AIC    LRT Pr(Chi)
<none>              38.66 455.63
Infl:Type:Sat 12    16.11 457.08  22.56 0.03175 *
Infl:Cont:Sat  4    37.47 462.44   1.19 0.87973
Type:Cont:Sat  6    28.26 457.23  10.41 0.10855
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>

>  pout1 <- polr(Sat ~ Infl + Type + Cont,
+      weights = Freq, data = housing)
>  summary(pout1)

Re-fitting to get Hessian

Call:
polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)

Coefficients:
                   Value Std. Error   t value
InflMedium     0.5663922 0.10465276  5.412109
InflHigh       1.2888137 0.12715609 10.135682
TypeApartment -0.5723552 0.11923800 -4.800107
TypeAtrium    -0.3661907 0.15517331 -2.359882
TypeTerrace   -1.0910073 0.15148595 -7.202036
ContHigh       0.3602803 0.09553577  3.771156

Intercepts:
            Value   Std. Error t value
Low|Medium  -0.4961  0.1248    -3.9740
Medium|High  0.6907  0.1255     5.5049

Residual Deviance: 3479.149
AIC: 3495.149
>

>  foo <- matrix(out1$fitted.values, ncol = 3, byrow = TRUE)
>  qux <- matrix(housing$Freq, ncol = 3, byrow = TRUE)
> bar <- apply(foo, 1, sum)
>  foo <- sweep(foo, 1, bar, "/")
>  bar <- pout1$fitted.values[seq(1, nrow(pout1$fitted.values)) %% 3 == 0, ]
> foo / bar
         Low    Medium      High
3  1.0452380 0.9041704 1.0312918
6  1.0132938 0.9751822 1.0072312
9  1.0473452 0.9115444 1.0183943
12 1.0456862 0.8861349 1.0269885
15 1.0387415 0.9118421 1.0319886
18 1.1130179 0.7982834 1.0541219
21 0.9184340 1.1729125 0.9638071
24 0.8897817 1.2058649 0.9406340
27 0.9571953 1.0990954 0.9709930
30 1.0012790 1.0305020 0.9495238
33 1.0011230 1.0141581 0.9811639
36 1.1035380 0.8391064 1.0310752
39 1.0006360 0.9915024 1.0053137
42 0.9512385 1.1022984 0.9717152
45 0.9637437 1.0741329 0.9875714
48 1.0188330 0.9464716 1.0243012
51 0.9937709 0.9987641 1.0053171
54 1.0442592 0.9140069 1.0208639
57 0.8738083 1.2414589 0.9353196
60 0.8363502 1.3209097 0.9018316
63 0.8885786 1.2678440 0.9369228
66 0.9799070 1.0721476 0.9659830
69 0.9679769 1.0784904 0.9701787
72 1.0562745 0.9358045 1.0066221
> dev <- 2 * sum(qux * log(foo / bar))
> dev
[1] 9.065433