[R] p values of plor
Prof Brian Ripley
ripley at stats.ox.ac.uk
Wed May 29 16:09:41 CEST 2013
AIC is a different story. To do hypothesis tests on terms, use anova()
or dropterm() (as done in the book):
library(MASS)
example(polr)
dropterm(house.plr, test = "Chisq")
Single term deletions
Model:
Sat ~ Infl + Type + Cont
Df AIC LRT Pr(Chi)
<none> 3495.1
Infl 2 3599.4 108.239 < 2.2e-16
Type 3 3545.1 55.910 4.391e-12
Cont 1 3507.5 14.306 0.0001554
That is an object of class "anova", so can be saved and subsetted.
On 28/05/2013 22:30, David Winsemius wrote:
>
> On May 27, 2013, at 11:05 PM, Prof Brian Ripley wrote:
>
>> On 28/05/2013 06:54, David Winsemius wrote:
>>>
>>> On May 27, 2013, at 7:59 PM, meng wrote:
>>>
>>>> Hi all:
>>>> As to the polr {MASS} function, how to find out p values of every
>>>> parameter?
>>>>
>>>>
>>>>> From the example of R help:
>>>> house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data =
>>>> housing)
>>>> summary(house.plr)
>>>>
>>>>
>>>> How to find out the p values of house.plr?
>>>
>>> Getting p-values from t-statistics should be fairly straight-forward:
>>>
>>> summary(house.plr)$coefficients
>>
>> And what distribution are you going to use to compute the p-values?
>
> I should have responded with my first impulse: "If the authors didn't provide p-values, then perhaps they don't think they are credible."
>
>>
>> Hint: there is no exact distribution theory for POLR fits and the asymptotic theory can be far enough off to be seriously misleading (just as for the two-class case, logistic regression: see MASS the book). That is why likelihood-ratio tests are recommended in MASS, not Wald tests.
>
> And so the more correct answer would be to use stepAIC? I would have thought sequential removal of terms with comparisons of deviance estimates might be informative. This is what I get with that data:
>
>> house.AIC.1 <- stepAIC(house.plr, list(upper=~., lower=~1) )
> Start: AIC=3495.15
> Sat ~ Infl + Type + Cont
>
> Df AIC
> <none> 3495.1
> - Cont 1 3507.5
> - Type 3 3545.1
> - Infl 2 3599.4
>>
>
> So something along those lines seems to be happening, but I am not able to extract those values programmatically, nor am I able to see how they even get displayed.
>
>> class(house.AIC.1)
> [1] "polr"
>> str(house.AIC.1$anova)
> Classes ‘Anova’ and 'data.frame': 1 obs. of 6 variables:
> $ Step : Factor w/ 1 level "": 1
> $ Df : num NA
> $ Deviance : num NA
> $ Resid. Df : num 1673
> $ Resid. Dev: num 3479
> $ AIC : num 3495
>
>
> Which lead me to look at:
>
> getAnywhere(print.polr)
>
> But that was uninformative to my level of reading R code. The AIC trials seem to get printed by stepAIC() but are not saved in the returned object.
>
> ---
>
> David Winsemius
> Alameda, CA, USA
>
--
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
More information about the R-help
mailing list