[R] deviance in polr method

[Prof Brian Ripley]

On Tue, 13 Jan 2009, David Winsemius wrote:

> I remember have the same consternation using GLIM with binomial models on 
> grouped and ungrouped data, but I was counseled by my betters only to 
> consider differences in models. The differences in deviance are the same up 
> to rounding error.
>
>> 859.8018 - 711.3479
> [1] 148.4539
>
>> 168.8-20.3
> [1] 148.5
>
> So I would ask if these really are different answers?

I think Gerard already said that, and accepts that the test between 
two models is the same.  He seems to be asking about using the 
comparison with one or other saturated model as some sort of test of 
model fit, something that is not normally recommended, including in 
the book for which polr() is support software.  In any case, not an R 
question any more.

>
> -- 
> David Winsemius
> Heritage Labs
>
> On Jan 13, 2009, at 10:06 AM, Prof Brian Ripley wrote:
>
>> Withut looking up your reference, are you not comparing grouped and 
>> ungrouped deviances?  And polr() does not say anything about accepting a 
>> model (or not), only about the comparison between two models.
>> 
>> 'Deviances' are in comparison with some 'saturated' model, and I would say 
>> that M&N are comparing with a separate fit to each group, polr() with 
>> correctly predicting each observation.  Both are valid, but refer to 
>> different questions.
>> 
>> On Tue, 13 Jan 2009, Gerard M. Keogh wrote:
>> 
>>> 
>>> Dear all,
>>> 
>>> 
>>> I've replicated the cheese tasting example on p175 of GLM's by McCullagh
>>> and Nelder. This is a 4 treatment (rows) by 9 ordinal response (cols)
>>> table.
>>> 
>>> Here's my simple code:
>>>
>>>    #### cheese
>>>    library(MASS)
>>>    options(contrasts = c("contr.treatment", "contr.poly"))
>>>    y = c(0,0, 1, 7, 8,8,19, 8,1, 6,9,12,11, 7,6, 1, 0,0, 1,1, 6, 8,23,7,
>>>    5, 1,0, 0,0, 0, 1, 3,7,14,16,11)
>>>    p =
>>>    c(1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9)
>>>    x1 =
>>>    c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
>>>    x2 =
>>>    c(0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
>>>    x3 =
>>>    c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0)
>>>    # catgeory 4 is used as baseline
>>>    ord.rp = ordered(p)
>>>    fit0.polr = polr(ord.rp ~ 1, weights=y)
>>>    fit1.polr = polr(ord.rp ~ x1 + x2+x3, weights=y)
>>>    summary(fit1.polr)
>>>    anova(fit0.polr,fit1.polr)
>>> 
>>> This works and "summary" gives the correct parameter estimates but I have 
>>> a
>>> problem with the deviance.
>>> Anova gives
>>>    Likelihood ratio tests of ordinal regression models
>>>
>>>    Response: ord.rp
>>>             Model Resid. df Resid. Dev   Test    Df LR stat. Pr(Chi)
>>>    1            1       200   859.8018
>>>    2 x1 + x2 + x3       197   711.3479 1 vs 2     3 148.4539       0
>>> 
>>> 
>>> In McCullagh the deviance for the null model is given as 168.8 on 24 d.f.
>>> and the fitted model is 20.3 on 21 d.f.
>>> This means the change in POLR and in the book is 148.5 on 3 d.f. so the
>>> model gives a significant improvement.
>>> 
>>> But the model deviance from POLR is 711 vs. 20.3 from McCullagh.
>>> Thus POLR says we should reject the model fit with chi-sq = 35 while
>>> McCullagh say we should accept it
>>> 
>>> Could someone explain what I should do when it comes to accepting or
>>> rejecting the model fit itself in R?
>>> 
>>> I'm using R 2.8.0.
>>> 
>>> Thanks.
>>> 
>>> Gerard


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595