[R] Strange parametrization in polr

BXC (Bendix Carstensen) bxc at steno.dk
Thu Jan 8 15:50:12 CET 2004


In Venables \& Ripley 3rd edition (p. 231) the proportional odds model
is described as:

logit(p<=k) = zeta_k + eta

but polr apparently thinks there is a minus in front of eta,
as is apprent below.

Is this a bug og a feature I have overlooked?


Here is the naked code for reproduction, below the results.
------------------------------------------------------------------------
---
version
library( MASS )
data( housing )
hnames <- lapply( housing[,-5], levels )
house.plr <- polr( Sat ~ Infl + Type + Cont, data=housing, weights=Freq
)
summary( house.plr )
newdat <- expand.grid( hnames[-1] )[1:5,]
cbind( newdat, predict( house.plr, newdat, type="probs" ) )
# Baseline probs:
diff( c(0,tigol( c(-0.4961,0.6907) ), 1) )
# First level of Infl:
diff( c(0,tigol( c(-0.4961,0.6907) + 0.5663922 ), 1) )
# But the change of sign for eta is needed to reproduce the fitted
values:
# Line 2:
diff( c(0,tigol( c(-0.4961,0.6907) - 0.5663922 ), 1) )
# Line 5:
diff( c(0,tigol( c(-0.4961,0.6907) - 0.5663922 + 0.5723552 ), 1) )
------------------------------------------------------------------------
---

Here is the resulting output:
------------------------------------------------------------------------
---
> version
         _              
platform i386-pc-mingw32
arch     i386           
os       mingw32        
system   i386, mingw32  
status                  
major    1              
minor    8.1            
year     2003           
month    11             
day      21             
language R              
> library( MASS )
> data( housing )
> hnames <- lapply( housing[,-5], levels )
> house.plr <- polr( Sat ~ Infl + Type + Cont, data=housing,
weights=Freq )
> summary( house.plr )

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.3661908 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 
> newdat <- expand.grid( hnames[-1] )[1:5,]
> cbind( newdat, predict( house.plr, newdat, type="probs" ) )
    Infl      Type Cont       Low    Medium      High
1    Low     Tower  Low 0.3784485 0.2876755 0.3338760
2 Medium     Tower  Low 0.2568261 0.2742125 0.4689614
3   High     Tower  Low 0.1436927 0.2110841 0.6452232
4    Low Apartment  Low 0.5190450 0.2605077 0.2204473
5 Medium Apartment  Low 0.3798522 0.2875967 0.3325511
> # Baseline probs:
> diff( c(0,tigol( c(-0.4961,0.6907) ), 1) )
[1] 0.3784576 0.2876650 0.3338774
> # First level of Infl:
> diff( c(0,tigol( c(-0.4961,0.6907) + 0.5663922 ), 1) )
[1] 0.5175658 0.2609593 0.2214749
> # But the change of sign for eta is needed to reproduce the fitted
values:
> # Line 2:
> diff( c(0,tigol( c(-0.4961,0.6907) - 0.5663922 ), 1) )
[1] 0.2568335 0.2742035 0.4689630
> # Line 5:
> diff( c(0,tigol( c(-0.4961,0.6907) - 0.5663922 + 0.5723552 ), 1) )
[1] 0.3798613 0.2875862 0.3325525
> 
------------------------------------------------------------------------
-----

----------------------
Bendix Carstensen
Senior Statistician
Steno Diabetes Center
Niels Steensens Vej 2
DK-2820 Gentofte
Denmark
tel: +45 44 43 87 38
mob: +45 30 75 87 38
fax: +45 44 43 07 06
bxc at steno.dk
www.biostat.ku.dk/~bxc




More information about the R-help mailing list